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1.  INTRODUCTION 


The  determination  of  the  gravity  field  is  one  of  the  major  tasks  of  geodesy  besides  precise 
positioning.  Especially  in  the  last  two  decades  great  progress  could  be  achieved.  On  the 
one  side  satellite  observations  became  available  which  contributed  significantly  to  the 
construction  and  improvement  of  global  earth  models  ( SZABO  1986).  The  best  of  those 
are  nowadays  complete  to  degree  and  order  m  =  n  =  360  (RAPP,  CRUZ  1986) 
corresponding  to  a  resolution  of  30’  *  30’  mean  gravity  anomalies,  say  50  *  50  km2. 
Satellite  altimetry  has  reversed  the  state  of  global  knowledge  of  the  gravity  field  with 
respect  to  land  and  oceans.  Whereas  still  in  the  1960’s  the  gravity  field  over  land  was 
much  better  known  than  over  sea,  now  the  opposite  became  true. 

On  the  other  side  integrated  data  processing  techniques  were  developed  which  are  able  to 
give  best  estimates  by  combining  all  available  geodetic  observations.  Their  field  of 
application  can  be  found  more  or  less  in  the  local  or  regional  range.  Least-squares 
collocation  and  integrated  geodesy  adjustment,  first  developed  and  proposed  by  KRARUP 
(1969),  EEG,  KRARUP  (1978),  are  the  key  developments  on  this  theoretical  and  data 
analysis  side. 

Despite  of  these  great  achievements  on  the  observation  recovery  and  the  analysis  side,  the 
demands  grow  further  for  a  better  resolution  and  accuracy  of  the  earth’s  gravity  field. 
With  respect  to  global  aspects  the  needs  are  coming  from  satellite  geodesy  where  most 
missions  require  an  orbit  determination  in  the  submeter-  or  even  decimeter  level.  They  are 
coming  from  oceanic  circulation  investigations,  sea-level  variability  studies  and  marine 
geophysics.  Crustal  dynamics  and  plate  tectonic  research  has  to  be  mentioned,  too.  On  the 
regional  or  local  scale  the  advent  of  the  Global  Positioning  System  has  started  a  revolution 
in  surveying  techniques.  Although  the  measurements  yield  threedimensional  positions,  the 
height  problem  cannot  be  solved  due  to  its  dynamical  character.  Thus,  precise  geoid 
heights  -  nothing  else  than  a  functional  of  the  geopotential  -  have  to  be  known  in  order  to 
derive,  for  example,  orthometric  heights  from  GPS  ellipsoidal  heights.  A  method  which 
can  lead  to  a  replacement  of  the  traditional  spirit  levelling,  so  far  physical  geodesy  is  able 
to  provide  the  geoid  with  sufficient  accuracy.  Prerequisite  again  is  the  precise  knowledge  of 
gravity  (anomalies)  and/or  vertical  deflections  in  those  areas. 

What  is  the  state-of-the  art  in  the  knowledge  of  the  gravity  field?  Although  due  to  satellite 
altimetry  nearly  a  complete  gravity  coverage  over  the  oceans  can  be  used  for  the 
construction  of  global  earth  models,  there  are  still  something  like  5500  1°  *  1°  blocks  over 
land  where  mean  gravity  values  for  certain  reasons  are  not  available.  This  corresponds  to 
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approximately  one  third  of  the  earth’s  land  surface  (see  Fig.  1.1).  Those  gaps  have  to  be 
filled  by  so-called  geophysical  anomalies,  that  means  that  the  gravity  anomaly  values  have 
to  be  estimated  on  the  basis  of  geophysical  information  and  models.  There  is  no  doubt  that 
an  improvement  in  those  procedures  and  gravity  values,  respectively,  results  in  more 
accurate  global  geopotential  models. 

Approximation  methods  like  collocation  are  able  to  combine  different  geodetic  and 
gravimetric  observations  to  determine  the  gravity  field  and  its  functionals  with 
wavelengths  smaller  than  the  resolution  of  the  global  earth  models.  Mainly,  gravity 
(anomalies)  and  astronomic  latitudes  and  longitudes  (deflections  of  the  vertical)  are 
combined.  For  gravity  anomaly  interpolation  a  certain  reduction  has  to  be  carried  out  in 
order  to  perform  the  real  numerical  computations  in  a  smooth  field,  often  referred  to  as 
"remove-restore''  technique.  For  the  determination  of  the  geoid  as  a  boundary  value 
problem,  gravity  and  other  data  have  to  be  reduced  down  to  the  geoid.  For  both  purposes 
mentioned  before  digital  terrain  models  have  to  be  available  and  rock  density  assumptions 
are  necessary.  Although  the  gravity  field  is  caused  by  mass  irregularities  within  the  earth  - 
see  Newton’s  law  of  gravitation  -  geodetic/gravimetric  algorithms  for  the  determination  of 
the  gravity  field  (prediction,  interpolation,  etc.)  do  not  use  geophysical  data,  like  density 
and  mass  distribution  to  construct  the  gravitational  potential, 

V(x)  =  G  fff  (1-1) 

n  ls  -  I' 

where  G  is  the  gravitational  constant,  ve  is  the  volume  of  the  earth  E,  p  is  the  density  and 
x,  y  are  threedimensional  position  vectors. 

The  imperfect  knowledge  of  the  density  distribution  within  the  earth,  however,  results  in 
the  fact,  that  the  gravitational  potential  outside  the  earth  has  to  be  determined  by  solving 
a  boundary  value  problem  with  gravity  anomalies  as  boundary  data,  for  example. 

For  both  methods,  smoothing  as  well  as  reduction,  density  assumptions  have  to  be  made. 
In  case  of  isostatic  gravity  anomalies  even  hypotheses  on  the  model  of  isostatic 
compensation  and  the  depths  of  the  Mohorovicic-discontinuity  are  necessary  to  perform 
the  necessary  computations.  Thus,  the  geophysical  (model)  data  are  used  indirectly  in 
geodesy.  With  increasing  observations  in  seismics,  geophysics  is  nowadays  able  to  give 
reasonable  density  values  or  seismic  velocity  information  for  whole  regions.  An  excellent 
example  is  the  digital  surface  density  model  of  WALACH  (198 7).  Often  also  density  values 
at  discrete  locations  are  available.  This  seems  to  be  now  the  time  to  think  about  the  direct 
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Fig  1.1.  Locations  of  the  6331  geophysically  predicted  1°  *  1°  mean  gravity  anomalies  A 
(File:  JUN86.TUG87.DMAMOE87,  RAPP  1989,  pers.  comm.). 
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use  of  geophysical  data  in  geodetic  approximation  methods  for  the  determination  of  the 
gravity  field.  First  proposals  in  the  last  years  are  reviewed  in  chap.  4.  Two  warnings  are 
always  coming  from  geophysicists  when  we  try  to  use  those  available  density  data  and 
establish  a  certain  mathematical/physical  relationship  between  density  and  gravity  (see 
also  chap.  2). 

Density  is  varying  with  depth  and  can  change  drastically  also  in  lateral  direction. 
Consequently,  it  is  very  hard  to  establish  a  density  model.  No  c*  ;bt,  it  has  to  be 
threedimensional.  The  need  of  geodesy,  to  generalize  the  density  information  -  both, 
laterally  in  the  form  of  certain  mean  blocks,  and  vertically,  in  the  form  of  a  minimum  of 
layers  -  implies  a  certain  loss  on  geophysical  information.  It  is,  however,  required  for  the 
sake  of  numerical  handling  of  the  problem  in  the  computer.  Insofar,  a  digital  model  of 
(surface)  rock  densities  has  its  problem.  What  does  surface  density  mean?  Obviously,  the 
first  some  hundred  meters  down  to  sea-level  are  characterized:  a  small  distance  compared 
to  the  earth’s  radius. 

Due  to  the  complexity  of  the  earth’s  crust  a  simple  mathematical  relationship  between 
density  p  and  gravity  with  sufficient  accuracy  exists  only  locally,  g  =  f(p).  However,  the 
problem  in  opposite  direction,  p  =  f(g)  is  ambiguous  and  usually  called  " inverse  problem " 
in  geophysics. 

Nevertheless,  (more  details  are  outlined  in  the  next  chapter  on  the  definition  of  the  prob¬ 
lem),  in  view  of  all  those  generalizations  and  approximations  which  have  to  be  made,  it  is 
better  to  use  such  density  information,  for  example,  in  gravity  (anomaly)  interpolation 
than  to  leave  the  gravity  gaps  white  or  to  use  globally  averaged  estimates  like  p  =  2.67  g 
c.m-3.  Even  probabilistic  approaches  are  justified  ("When  in  doubt,  smooth",  Sir  Harold 
Jeffreys). 

This  report  proposes  some  first  (simple)  ideas,  proposals,  and  approaches  how  to  use  geo¬ 
physical  data  in  the  determination  of  the  earth’s  gravity  field  and  its  functionals  (gravity, 
deflections  of  the  vertical,  geoid  heights).  Numerical  investigations  with  real  data  are  espe¬ 
cially  carried  out  in  the  field  of  gravity  interpolation. 
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2.  DEFINITION  OF  THE  PROBLEM 


We  will  assume  that  we  have  made  the  following  measurements  on  the  earth’s  surface  (or 
in  space), 

-  gravity  g,  astronomic  latitude  $,  astronomic  longitude  A,  gravity  gradients 
d2W/3xidxj  (second-order  derivatives  of  the  gravity  potential)  or,  in  other  form,  after 
subtracting  reference-  or  model-{ normal)  values: 

gravity  anomalies  Ag,  components  of  deflections  of  the  vertical  in  north-south  and  in 
east-west  direction  (£,t?),  resp.,  anomalous  gravity  gradients  ftT/dxidxj  . 

We  further  assume,  that  geophysical  information  is  available,  in  detail, 

-  density  data,  characterizing  as  discrete  point  or  mean  block  values  over  a  certain 
lateral  area  a  specified  layer  of  the  earth’s  crust, 

-  seismic  velocities  v,  referring  as  representative  values  to  certain  portions  or  layers 
within  the  earth’s  crust, 

-  seismic  displacements. 

In  the  ideal  case  density  and  seismic  velocities  are  given  in  the  form  of  a  threedimensional 
digital  model. 

We  are  looking  for  a  solution  for  the  gravity  potential  W  (or  disturbing  potential  T) 
outside  the  earth  by  combining  all  available  geodetic  and  geophysical  information  as 
mentioned  above.  From  this  representation  of  W  (or  T)  we  like  to  derive  all  functionals, 
like  gravity  anomalies,  deflections  of  the  vertical,  geoidal  heights,  etc. 

A  good  example  for  the  problem  stated  above  in  the  real  world  would  be  the  prediction  of 
gravity  anomalies  in  unsurveyed  areas  (without  gravity  observations)  where,  however, 
density  estimates  and/or  seismic  velocities  either  from  geophysical  investigations  are 
available,  or  where  a  reasonable  guess  of  the  density  variation  based  on  geological  features 
can  be  made.  In  addition,  gravity  values  in  the  surrounding  areas  are  available  which 
could  be  considered  in  a  combination  solution. 

Central  point  in  such  a  solution  algorithm  is  the  mathematical-physical  foundation  of  the 
relationship 

seismic  velocity  v  •— >  density  p  •— *  gravity  g  (or  gravity  potential  W)  . 
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It  was  already  mentioned  that  the  determination  of  internal  masses  from  (surface)  gravity 
data,  p  =  f(g),  has  no  unique  solution  and  is  treated  in  geophysics  as  so-called  inverse 
problem. 


<=> 


GEOPHYSICS 


GRAVITY  FIELD  OF  THE  EARTH 

It  i 

DENSITY 

P 

It  T 

SEISMICS 

seismic  velocities 


Fig.  2.1.  Gravity,  density,  and  seismic  data  in  geodesy  and  geophysics  (the  corresponding 
arrows  indicate  the  main  goals  and  working  directions  of  geodesy  (=>)  and  geo¬ 
physics  (->)). 

There  are,  however,  two  major  differences  to  geophysics  in  our  task  which  are  also 
illustrated  in  Fig.  2.1.  First,  we  are  only  interested  in  the  solution  in  the  one  direction  v  — > 
p  — ♦  g.  A  fully  established,  mathematically  and  physically  correct  inverse  relationship  or 
solution  is  not  of  primary  interest  for  us,  since  we  are  only  interested  in  the  improvement 
of  the  exterior  gravity  field.  It  may  be  possible,  for  example,  that  a  certain  relationship 
between  p  and  g  leads  to  a  reasonable  solution  for  g  =  f (/?);  however,  the  inverse-relation 
applied  in  the  algorithm  may  yield  unacceptable  results.  The  second  point  which  should  be 
raised  here  when  defining  the  problem  is  the  following.  If  we  are  faced  with  the  problem 
that  no  gravity  observations  are  available  in  certain  regions  of  the  world  even  a  guess  of 
the  geological/geophysical  situation  may  help  in  the  prediction  of  smoothed  gravity  values 
and  may  be  better  than  to  use  zero-values  in  the  spherical  harmonics  solution  of  an  earth 
model.  It  is  obvious,  that  for  such  a  task  some  kind  of  generalization  has  to  take  place 
which  may  drastically  smooth  local  geophysical  features. 
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Finally,  we  have  to  decide  on  the  character  of  the  relationship  v  — »  p  — ►  g.  Should  we 
define  it  in  the  deterministic  way  or  is  it  opportune  to  apply  the  concept  of  random 
processes  in  such  a  (prediction)  algorithm.  There  are  good  reasons  to  vote  for  the  last. 
Gravity  field  determination  involves  functions,  or  better  functionals  like  geoid  heights, 
gravity  disturbances,  deflections  of  the  verticals,  etc.  which  have  to  be  evaluated  from  the 
same  source.  Our  knowledge  about  the  interior  of  the  earth  and  the  needed  mathematical/ 
physical  relationships  are  limited.  Thus,  the  concept  of  probability  enters  to  some  extent 
into  our  problem.  Even  from  the  (visual)  distribution  of  the  masses  within  the  earth  one 
might  think  of  some  irregular  distribution  which  justifies  the  application  of  stochastic  or 
random  theory.  No  doubt,  that  it  would  also  extend  the  already  applied  concept  of 
geodetic  collocation  in  integrated  geodesy.  These  arguments  all  have  led  us  to  the  decision 
to  look  -  besides  other  deterministic  methods  -  also  for  an  (overall)  stochastic  concept. 
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3.  THE  INFLUENCE  OF  SEISMIC  STRUCTURES  ON  THE  GRAVITY  FIELD: 
THE  GEOPHYSICIST’S  VIEW 


In  order  to  use  geophysical  data  for  gravity  interpolation  we  have  to  understand  more  in  de¬ 
tail  the  underlying  mathematical/physical  relationship  between  the  corresponding  observa¬ 
tions.  Since  those  relations  often  have  model  character  or  represent  simplified  physics  it  is 
worthwhile  to  verify  them  also  with  real  data.  Therefore  this  chapter  mainly  tries  to  inves¬ 
tigate  available  models  and  includes  also  some  inverse  computations  determining  density 
contrasts  between  the  different  layers  of  the  earth’s  crust. 

The  layer  boundaries  determined  from  seismic  or  seismological  data  are  the  most  direct 
evidence  for  physical  inhomogeneity  within  the  earth’s  crust  and  upper  mantle.  Since 
density  tends  to  correlate  with  seismic  velocity,  undulating  layer  boundaries  are  very  likely 
to  influence  the  gravity  field.  The  recently  released  digital  density  models,  see,  e.g.,  WA- 
LACH  (1987),  represent,  however  only  surface  densities.  In  spite  of  their  rather  exact  know¬ 
ledge  they  cannot  be  reliably  continued  downward.  More  knowledge  about  the  density  dis¬ 
tribution  at  depth  is  necessary,  if  we  want  to  compute  and  improve  the  gravity  field  and  its 
functionals. 

The  most  important  and  widely  accepted  way  to  determine  the  rock  densities  at  depth  is  to 
apply  the  velocity-density  systematics  or  relationships.  They  are  based  on  experimental 
laboratory  measurements  (e.g.,  BIRCH  I960,  1961;  LIEBERMANN,  RINGWOOD  1973)  or 
on  solid  state  physical  theory  (e.g.,  BIRCH  1979;  ANDERSON  1973).  Another  approach  is 
the  inversion  of  gravity  data  with  the  assumption  of  "known"  geometrical  structure,  e.g., 
from  seismic  studies,  and  the  additional  assumption  that  bodies  of  constant  seismic  velocity 
are  also  uniform  in  density  (e.g.,  JACOBY  1973,  1975). 

We  must  note,  however,  that  gravity,  as  any  potential  field,  cannot  be  inverted 
unambiguously  but  that  in  principle  there  is  an  infinite  number  of  solutions.  Nevertheless 
the  number  can  be  drastically  reduced  by  adding  information,  e.g.,  from  geology,  and 
seismic  or  other  geophysical  information.  If  the  structure  of  the  crust  is  sufficiently  known, 
the  gravity  effect  of  finite  parts,  e.g.,  bodies  or  layers,  can  be  computed  at  any  surface 
points.  If  only  their  densities  are  unknown,  one  may  determine  them  with  the  aid  of  the 
least-squares  method  by  fitting  the  computed  sum  of  the  individual  bodies  to  the 
observations.  Subsequently  one  may  relate  the  density  variation  found  with  the  already 
known  seismic  velocities  to  determine  the  velocity-density  systematics.  From  this,  it  is 
hoped,  the  physical  state  and  the  chemical  and  petrological  composition  of  the  earth’s 
interior  will  be  better  understood. 
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In  principle,  it  should  be  possible  to  extend  this  method  to  gravity  prediction.  If  the  gravity 
field  is  known  incompletely,  e.g.,  only  in  parts  of  a  region  or  only  in  its  long  wavelength 
components,  but  published  seismic  models  are  fairly  detailed,  one  may  attempt  the 
following  approach.  First  find  the  best  fitting  densities  as  sketched  above  and  then  compute 
the  gravity  anomalies  and/or  geoid  undulations  that  are  caused  by  the  seismic  structure 
with  the  density  variation  determined.  However,  whether  this  approach  is  successful  has  to 
be  tested  with  practical  examples. 


3.1  Computation  of  the  gravity  effect  of  given  model  bodies 

The  starting  point  for  the  computation  of  the  gravity  effect  of  given  masses  is  the 
contribution  of  an  infinitesimal  mass  element  dm  to  the  gravitational  potential  0V .  It  can 
be  expressed  as: 

w(5)  =  G  fff  Ts^TT  (3~1) 

V 

where  G  is  the  universal  gravitational  constant  ((6.673  ±  0.003  )  *  10‘8  g"1  cm3  s-2),  and  the 
vectors  y  =  {^,r/,(}  and  x  =  {x,y,z}  represent  the  coordinates  of  the  mass  element  and 
the  observation  point  (where  the  gravity  effect  is  to  be  calculated),  respectively  (see  Fig. 
3.1). 

Differentiating  equation  (3-1)  with  respect  to  z  (vertical  downward)  assuming  a  local 
cartesian  threedimensional  coordinate  system  and  expressing  |x-y|  by  their  components 
renders  the  gravity  effect  fig  of  the  mass  body  (density  is  denoted  by  p  ) 

«*)  -  “  G  fJJ,  (3'2) 


The  integral  (3-2)  can  be  most  easily  computed  for  two-dimensional  masses.  Such  masses 
can  be  represented  by  a  polygon  (with  j  corners)  in  the  normal  section  of  the  mass  which  is 
assumed  to  be  infinite  horizontally  at  the  right  angle  to  the  section.  According  to  JUNG 
(1961)  the  gravity  effect  at  the  point  P  of  a  triangular  section  with  apex  P(x)  is  (see  Fig. 
3.2  for  the  notation  and  corresponding  quantities): 

fig{x)  =  2  Gp  (n  -  x)  sina  Jsina  In  gin  [a  +  C0SQ  “  WO]  •  (3“3) 
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The  effect  of  a  two-dimensional  mass  of  an  arbitrary  polygonal  cross  section  is  calculated 
by  dividing  the  polygon  into  j  triangles  between  two  neighbouring  corner  points  and  P  as 
the  apex.  The  total  gravity  effect  is  the  sum  of  all  triangle  effects  where  the  sign  changes 
along  the  circumference  of  the  polygon  such  that  the  effect  of  "external"  triangles  are 
subtracted  and  that  the  "interior"  area  (volume)  is  left  (see  Fig.  3.3). 

For  three-dimensional  masses,  more  appropriate  in  reality,  a  similar  approach  is  available. 
The  surface  of  a  body  is  divided  into  a  finite  number  of  triangles  or  one  approximates  the 
body  by  a  polyhedron.  Each  triangle  forms  a  tetrahedron  with  the  observation  point  P.  As 
in  the  two-dimensional  case,  the  gravity  effect  of  the  whole  body  is  the  sum  of  the  effects  of 
all  tetrahedra  with  a  similar  sign  rule.  Therefore  one  has  again  to  follow  strict  rules  when 
defining  the  orientation  of  the  triangles  at  the  body  surface,  e.g.,  by  surrounding  them 
always  clockwise  such  that  the  surface  normal  always  points  out  of  the  body.  The  analytical 
solution  of  the  gravity  integral  (3-2)  for  this  case,  i.e.,  the  gravity  effect  of  just  one 
tetrahedron  is  complicated  and  would  fill  several  pages  (QAV^AK  1988).  (A  computer  code 
in  FORTRAN77  is  available  to  the  authors.) 


3.2  Computation  of  geoid  undulations  caused  by  model  mass  anomalies 

For  the  computation  of  potential  anomalies  or  geoid  undulations  the  solution  has  been 
described  in  detail  by  CHAPMAN  (1979).  It  is  extendable  to  the  gravity  effect  without 
much  difficulty.  The  gravity  effect  is  obtained  by  differentiating  the  potential  solution  with 
respect  to  z.  Computing  the  geoid  means  to  determine  the  shape  of  an  equipotential  surface. 
Roughly  speaking  it  is  a  rotational  ellipsoid  with  small  perturbations  (undulations),  which 
are  related  to  the  disturbing  potential.  If  the  disturbing  potential  T  is  known,  the  geoid 
undulation  N  is  found  from  Bruns’  formula  (HEISKANEN,  MORITZ  1967): 

N=|,  (3-4) 

where  7  is  the  normal  gravity  value  at  the  reference  ellipsoid  vertically  below  (above)  the 
observation  point.  The  potential  5T  of  a  disturbing  point  mass  m  at  a  distance  1  is  given  by: 


(3-5) 
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If  this  definition  is  extended  to  the  expression  (3-1)  with  the  coordinates  used  there  (Fig. 
3.1)  the  integral  expression  for  the  geoid  becomes 


N(x.y,»)  =  g  fJJ,  if  *>  K 


(3-6) 


3.3  Estimation  of  the  normal  values  for  crustal  parameters 

In  gravimetric  modelling  of  the  crust  it  is  necessary  to  extend  bodies,  as,  e.g.,  layers  not 
limited  to  the  study  area,  beyond  its  margin  in  order  to  avoid  edge  effects.  One  may  use  a 
coarse  grid  to  represent  such  "external"  bodies  if  known,  taking  advantage  of  the  "remote 
zone  effect".  If  unknown,  as  often  the  case,  the  "external"  crustal  structures  have  to  be 
extended  anyway,  but  one  may  assume  standard  or  normal  models.  Deviations  of 
boundaries  from  their  normal  values  are  the  cause  of  gravity  and  geoid  anomalies  and  their 
knowledge  is  thus  important  as  a  reference.  For  gravity  calculations  the  crustal  model  must 
extend  several  hundred  kilometers  beyond  the  limits,  200  -  500  km  are  generally  sufficient. 
For  geoid  calculations  the  external  extension  must  be  greater,  at  least  600  -  1000  km, 
because  of  the  slower  decrease  of  the  potential  with  distance. 

In  estimating  the  normal  crustal  values,  e.g.,  of  the  Moho  depth,  one  may  make  use  of  a 
standard  functional  relationship  between  elevation  of  the  earth’s  surface  and  Moho  depth. 
This  is  because  the  isostatic  equilibrium  of  sufficiently  large  regions  of  the  earth’s  crust  is 
generally  expected  to  exist.  We  may  assume  the  average  linear  relationship,  which  has  been 
reported  by  TURCOTTE,  SCHUBERT  (1982),  based  on  the  Airy  model  of  isostasy  and  the 
assumption  of  standard  values  of  crustal  thickness  and  crust  as  well  as  mantle  densities. 

Z.  [km]  =  33.0  +  6.60  *  H  [km]  .  (3-7) 

Zm  is  the  depth  of  the  Mohorovidc  discontinuity  and  H  is  the  average  topographic  elevation 
of  a  compartment. 

If  one  has  sufficient  data  for  the  study  area,  one  may  determine  regionally  more  appropriate 
values  of  A  and  B  in  the  model  ZB  =  A  +  B  *  H  by  linear  regression  analysis.  In  order  to 
determine  this  isostatic  effect  and  the  normal  average  crustal  thickness  one  limits  the 
investigation  to  elevations  H  >  0.2  km.  Fig.  3.4  and  Table  3.1  present  the  results  of  such  an 
analysis  on  the  basis  of  all  compartments  in  the  study  area  except  those  in  the  southern 
Rhine  Graben,  where  the  Moho  is  anomalously  shallow,  and  the  southeast  part  of  the  study 
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area  where  our  Moho  data  are  merely  interpolated.  Table  3.1  also  includes  the  results  of 
GEISS  (1987)  for  comparison. 


Fig.  3.4.  The  correlation  of  Moho  depth  and  topographic  elevation  in  the  test  area 
"European  Alps". 


Table  3.1.  Regression  of  topography  and  crustal  thickness 
Zm  [km]  =  A  +  B*H  [km]. 


n  ...  number  of  data  dupletts, 

AA,AB  ...  standard  deviation  of  the  coefficients  A,B, 

R  ...  correlation  coefficient, 

STDV  ...  standard  deviation  of  the  regression  [km] 

Model  1  ...  crust  (H  >  0.2  km)  Central  Europe  and  Alps  (GE*SS  1987) 
Model  2  ...  Oceans  (H  <  -  0.2  km)  (GEISS  1987) 

Model  3  ...  Alps,  all  data  dupletts  (H  >  0.2  km) 

Model  4  ...  as  Model  3,  but  without  anomalous  and  uncertain  values 
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The  results  are  rather  similar,  but  the  correlation  coefficients  are  a  little  larger,  and  the 
standard  deviations  are  smaller  than  those  of  GEISS  (1987),  whose  investigation  was  based 
on  much  bigger  compartments  (12*  *  20’  and  1°  *  1°).  The  regression  shows  also  that  the 
relationship  between  elevation  and  crustal  thickness  fit  the  data  with  remarkably  little 
scatter  in  the  Alpine  region.  An  "isostatic"  Mohorovi£ic  discontinuity  computed  on  the 
basis  of  the  models  2  and  4  (Table  3.1)  shows  a  good  agreement  with  the  seismically 
determined  one  (Fig.  8.13),  in  fact,  it  agrees  within  the  reliability  limits  (Figs.  8.14  and 
8.15).  Note,  however,  that  the  dependence  of  Moho  depth  on  elevation  is  greater  in  the  seis¬ 
mic  model  than  in  the  theoretical  one  of  equation  (3-7),  based  on  Airy-isostasy  and  stand¬ 
ard  densities. 

Similar  regression  analyses  for  the  other  boundaries  (seismic  basement  Zb  and  lithosphere/ 
asthenosphere  boundary  Zi)  are  not  possible  or  sensible.  The  basement  does  not  show  any 
unique  correlation  with  the  topography  and  it  is  not  the  result  of  isostatic  equilibrium.  The 
same  is  true  for  the  lithospheric  root  which  may  be  a  relic  of  subduction.  Isostasy  is  the 
vertical  mass  equilibrium,  or  in  other  words,  the  tendency  towards  hydrostatic  pressure  at 
some  depth  with  a  laterally  varying  density  distribution  above  the  surface  of  pressure 
equilibrium.  Experience  with  other  large  scale  models  (JACOBY  1973,  1975)  have  shown 
that  this  surface  must  lie  below  the  deepest  Moho  depth,  but  it  is  not  identical  with  the 
lithosphere/asthenosphere  boundary.  Thus,  in  estimating  the  normal  Zb  and  Zj  values  one 
can  only  rely  on  the  literature. 

The  work  of  MOSTAANPOUR  ( 1 984 )  indicates  that  the  depth  of  the  seismic  basement  in 
large  parts  of  Europe  is  close  to  6.0  ±0.1  km.  The  mean  value  of  Zb  in  our  study  area  is 
only  5.5  km.  The  normal  thickness  of  the  "lid"  or  the  thickness  of  the  uppermost  mantle 
part  of  the  lithosphere  between  Moho  and  top  of  the  asthenosphere  is  given  by 
SUHADOLC,  PANZA  (1988)  for  the  Mediterranean  -  Alpine  region  to  be  70  -  80  km.  With 
the  normal  crustal  thickness  added  we  get  about  110  km  for  the  normal  lithosphere 
thickness.  For  our  purposes  this  value  is  questionable  since  we  are  interested  in  the 
anomalous  region  of  the  Alps.  The  lithosphere  thickness  in  marine  areas  (e.g.,  of  the 
Mediterranean)  is  only  40  -  60  km,  and  in  northern  central  Europe  the  values  are  also 
rather  low.  If  the  marginal  regions  of  our  models  are  assumed  to  follow  the  above  normal 
values,  the  error  introduced  into  the  gravity  computations  is  neglegibly  small  because  of  the 
rapid  gravity  decrease  with  distance.  In  calculating  the  geoid,  decreasing  as  1/1,  the  errors 
will  still  be  noticeable.  Generally,  a  given  mass  element  (of  the  upper  and  lower  crust  as 
well  as  the  uppermost  mantle  typical  for  an  Alpine  model,  each  with  the  same  lateral 
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extent,  directly  below  a  station)  at  13,  39  and  65  km  depth  has  effects  that  vary  as  1  :  0.3  . 
0.02  for  gravity  and  1  :  1  :  0.1  for  the  geoid;  for  horizontally  displaced  elements  the  geoid 
effect  decreases  even  less  with  depth.  Table  3.2  summarizes  the  mean  values,  reference 
depth  for  computed  values  and  depth  of  boundaries  used  as  extensions  beyond  the  model 
limits  for  the  following  gravimetric  computations. 


Table  3.2.  Mean  values  and  reference  depth  of  layer  boundaries. 


Layer 

MEAN 

REF 

C*H 

' _ 

Basement 

Mohorovicic  discontinuity 
Lithosphere  thickness 

5.53 

36.60 

143.19 

6.00 

30.00 

110.00 

MEAN  ...  mean  value  for  study  area  [km], 

REF  ...  reference  depth  [km], 

F  ...  extension  depth  outside  study  area  [km]. 


3.4  Two-dimensional  model  computations 

Although  two-dimensionality  is  a  gross  deviation  from  the  real  geometry  of  the  earth, 
two-dimensional  models  for  profiles  at  right  angles  to  the  main  strike  of  the  geological 
structures  are  a  simple  and  important  means  to  estimate  the  density  contrasts  at 
boundaries.  They  can  thus  give  an  indication  of  how  good  the  seismic  models  are.  The 
advantage  of  two-dimensional  modelling  is  that  it  requires  much  less  computing  time  than 
three-dimensional  models.  If  the  depth  extent  is  small  as  compared  to  the  approximately 
linear  horizontal  extent,  the  errors  caused  by  two-dimensionality  are  small.  In  the  study 
area  "European  Alps"  a  profile  was  chosen  that  is  perpendicular  to  the  main  axis  of  the 
Alps  (see  Fig.  3  5).  The  two-dimensionality  is  acceptable  in  view  of  the  maximum  model 
depth  of  less  than  200  km  and  the  radius  of  curvature  of  the  Alpine  axis  of  more  than  400 
km.  The  location  of  the  profile  is  determined  largely  by  the  availability  of  gravity  data  (see 
chap.  8.2.2)  and  parallels  the  Swiss  Geotraverse.  The  data  bases  are  the  area  means  of  the 
Bouguer  anomaly,  i.e,  a  smoothed  version  of  the  gravity  profile  (see  Fig.  3.5).  The  profile 
crosses  the  Ivrea  body,  which  is  possibly  a  high-density  piece  of  oceanic  crust  and  upper 
mantle  (lithosphere)  obducted  onto  continental  lithosphere  from  the  south.  The  Ivrea  body 
was  approximated  by  an  additional  model  body  dipping  to  the  south  as  indicated  in  the 
seismic  data  (BERCKHEMER  1968).  In  fitting  gravity,  this  assumption  proved  to  be  very 
important  to  reduce  the  residuals  or  "errors"  and  to  improve  the  density  estimation  for  the 
three  layers  of  our  general  model,  shown  in  Fig.  3.6  as  a  cross  section  through  the 
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three-dimensional  models  (Figs.  8.11,  8.13,  8.16).  Our  experience  is  that  the  neglect  of 
important  geometrical  model  aspects  will  adversely  influence  the  density  estimates  of  the 
rest  of  the  model. 

In  computing  the  densities  we  advanced  from  simple  models  to  as  complete  ones  as  possible. 
In  the  sequence,  we  first  (model  1)  asssumed  only  the  Mohorovicic  discontinuity  to  exist, 
then  we  added  the  seismic  basement  (model  2)  and  finally  completed  the  model  by  adding 
the  lithosphere/asthenosphere  boundary.  The  Ivrea  body  was  always  included.  The  results, 
i.e.,  our  estimates  of  the  various  density  contrasts  (across  the  boundaries  relative  to  the 
mantle  below  the  lithosphere)  are  presented  in  Table  3.3. 


Table  3.3.  Density  contrasts  from  two-dimensional  model  calculations  (fit  of  Bouguer 
anomaly)  relative  to  sublithospheric  mantle. 


Model 

Co 

mm 

c2 

c3 

C4 

R 

STDV 

1 

233.0±21.6 

-0.24*0.02 

0.13*0.05 

RSI 

2 

310.3*30.9 

-0.41*0.06 

-0.27*0.02 

0.12*0.04 

0.95 

KS6I 

3 

129.2*23.5 

-0.40*0.03 

-0.27*0.01 

0.06*0.02 

0.09*0.01 

0.99 

7.45 

C0  ...  constant  gravity  value  [mGal],  only  formal, 

C!  ...  density  contrast  of  lower  crust  [g/cm3], 

C2  ...  density  contrast  of  upper  (sedimentary)  crust  [g/cm3], 

C3  ...  density  contrast  of  lithosphere  fg/cm31, 

C4  ...  density  contrast  of  Ivrea  body  [g/cm3], 

R  ...  correlation  coefficient 

STDV  ...  standard  deviation  of  model  gravity  anomaly  [mGal], 

Model  1  ...  Moho  and  Ivrea  body  included, 

Model  2  ...  Moho,  seismic  basement,  Ivrea  body, 

Model  3  ...  Moho,  seismic  basement,  Ivrea  body,  and  lithosphere  /  as- 
thenosphere  boundary. 


The  results  demonstrate  that  the  models  get  better  if  more  details  (of  seismic  origin)  are 
included,  particularly  the  lithosphere/asthenosphere  boundary;  the  fit  is  greatly  improved 
(standard  deviation  less  than  half  from  model  1  to  3),  and  the  density  contrasts  assume 
values  that  are  very  plausible  in  view  of  seismic  velocities.  If  a  reference  density  of 
3  3  g/cm3  (KISSLING  1982)  is  assumed  for  the  uppermost  part  of  the  earth’s  mantle  below 
the  Alps,  the  densities  are  similar  to  the  densities  in  a  north-south  profile  (part  of  the 
European  Geotraverse)  which  was  computed  by  SCHWENDENER,  MULLER  (1985).  They 
found  the  following  values:  upper  crust  assumed  and  fixed  at  2.73  g/cm3;  middle  crust 
2.86  g/cm3,  lower  crust  2.957  g/cm3,  and  uppermost  mantle  3.264  g/cm3.  In  our  own 
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Fig  3. 


Fig.  3. 


Fig.  3.7. 


The  investigated  two-dimensional  profile  parallel  to  the  Swiss  geotraverse. 


Seismic  cross  section  along  the  two-dimensional  profile,  (1)  upper  crust 
above  seismic  basement,  (2)  lower  crust  between  seismic  basement  and 
Moho,  (3)  lithosphere,  (4)  Ivrea  body 


Observed  and  calculated  gravity  along  the  two-dimensional  profile  in 
[raGalj. 
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two-layer  model  the  corresponding  values  are:  upper  crust:  2.845  g/cm3,  middle  and  lower 
crust:  2.973  g/cm3,  and  uppermost  mantle  3.30  g/cm3  (assumed) 

The  positive  density  contrast  from  the  asthenosphere  to  the  lithosphere  is  expected  from 
geodynamic  arguments  and  from  the  low  seismic  velocities  in  the  asthenosphere.  Similar 
values  have  been  found  by  JACOBY  (1973,  1975)  and  by  SUHADOLC,  PANZA  (1988)  for 
the  Alps,  who  investigated  gravimetically  the  sub-Alpine  lithospheric  roots  and  obtained 
+  0.05  g/cm3.  In  contrast  GEISS  (1987)  found  a  value  of  about  -  0.01  g/cm3  which  is 
insignificant,  but  surprising. 

The  density  contrast  of  -  0.327  g/cm3  from  the  mantle  to  the  crust  falls  into  the  frequently 
assumed  range  of  -  0.30  to  -  0.35  g/cm3  (KISSLING  1982).  For  the  Ivrea  body  we  found 
the  density  to  be  nearly  that  of  the  upper  mantle.  This  was  the  result  of  a  free  fit  of  model 
gravity  to  the  observations.  These  results  give  us  some  confidence  into  the  model  and  we 
believe  that  this  is  an  argument  for  the  high  quality  of  seismic  data.  There  are  still  relative 
large  residuals  in  the  region  of  the  Ivrea  body  and  the  northern  Molasse  (see  Fig.  3.7);  these 
are  probably  caused  by  shallow  mass  anomalies  not  sufficiently  well  incorporated  into  our 
model.  The  deviations  from  the  two-dimensionality  assumed  will  also  contribute  to  the 
residuals. 


3.5  Three-dimensional  model  computations 

Two-dimensional  models  are  always  rather  arbitrary  and  incomplete,  though  very  simple. 
The  above  results  are  therefore  to  be  verified  by  computing  a  three-dimensional  model  for 
the  whole  study  area  "European  Alps".  We  computed  both  gravity  effects  and  geoid 
undulations  on  the  basis  of  the  shapes  of  the  boundaries  Z b,  Zm,  and  Z\.  In  computing  the 
geoid  we  must  also  take  into  account  the  effect  of  the  topography.  For  simplicity  we 
assumed  the  density  of  the  "topographic  masses"  to  be  2.67  g/cm3.  The  data  to  be  fitted 
were  the  geoid  undulations  represented  by  a  global  spherical  harmonic  expansion  to  degree 
and  order  360  (spatial  resolution  about  0.5  degree).  For  the  gravity  fit  we  again  used  the 
area  means  of  the  Bouguer  anomaly  (see  Chapter  8.2.2).  Compartments  near  the  Ivrea  body 
were  excluded  from  density  estimation  by  least-squares  gravity  fit  so  that  we  could  avoid 
the  introduction  of  an  additional  complicated  and  uncertain  three-dimensional  model  of  this 
structure.  The  geoid  representation  also  includes  no  effect  of  the  Ivrea  body. 

The  results  of  gravity  fitting  on  the  basis  of  1660  compartments  are  presented  in  Table  3.4. 
The  standard  deviations  are  much  greater  than  in  the  2-d  cases  The  density  contrasts 


18 


Table  3.4.  Density  contrasts  at  the  layer  boundaries  from  three-dimensional  model 
computations  (fitting  of  gravity). 


Model 


-28.23  ±1.39 
-15.96  ±1.48 
-36.39  ±4.80 


C, 


-0.191  ±0.011 
-0.176  ±0.012 


-0.207  ±0.005 
-0.253  ±0.005 
-0.278  ±0.008 


0.029  ±0.006 


& 


0.71 

0.76 

0.77 


STDV 


33.94 

31.37 

31.19 


Co  ...  constant  of  fitting  [mGal], 

Ci  ...  density  contrast  at  the  seismic  basement  [g/cm3], 

C2  density  contrast  at  the  Moho  fg/cm3], 

C3  ...  density  contrast  at  the  lithospnere  /  asthenosphere  boundary 

[g/cm3], 

R  ...  correlation  coefficient 

STDV'  ...  standard  deviation  of  the  model  [mGal], 

Model  1  ...  only  Moho  included, 

Model  2  ...  Moho  and  the  seismic  basement  included, 

Model  3  ...  Moho,  basement,  and  lithosphere  /  asthenosphere  boundary. 


computed  are,  however,  much  the  same.  The  density  contrasts  relative  to  the 
sub-lithospheric  mantle  and  the  absolute  values  (if  3.3  g/cm3  is  assumed  for  the 
lithosphere)  are:  lithosphere  +  0.029  (3.30)  g/cm3;  lower  crust  above  the  Moho:  -  0.249 
(3.022)  g/cm3;  and  upper  crust  above  the  seismic  basement:  -0.425  (2.846)  g/cm3.  Note 
that  the  density  contrast  between  the  lithosphere  and  asthenosphere  is  much  smaller,  only 
about  one  half  of  that  computed  for  the  two-dimensional  profile.  The  introduction  of  the  li- 
thosphere/asthenosphere  boundary  does  also  no  longer  reduce  the  standard  deviation  that 
much.  It  is  suspected  that  the  model  assumptions  -  mainly  the  homogeneity  of  the  layers 
between  the  boundaries  -  are  not  exactly  true  and  that  this  becomes  more  serious  the  larger 
the  region  included  in  the  modelling  (2-d  proGle  towards  3-d  model).  This  may  be  especially 
so  in  the  case  of  the  lithosphere  and  asthenosphere,  perhaps  the  lithospheric  roots.  The  two- 
dimensional  profile  exactly  crosses  the  root  below  the  western  Alps.  In  other  regions  the 
density  contrast  may  be  lower;  an  indication  is  also  the  result  of  GEISS  (1987)  -  see  above. 

The  effects  of  the  three  boundaries  with  the  density  contrasts  of  model  3  and  the  computed 
effect  resulting  from  the  fit  of  the  original  data  are  shown  in  Figs.  3.8  a,b,c  and  3.9.  The 
computed  gravity  map  for  the  region  of  Switzerland  is  already  quite  similar  to  that  of 
KISSLING  (1982)  in  which  the  effects  of  the  Ivrea  body  and  of  the  superficial  sediments  at 
the  northern  Alpine  margin  have  been  removed  from  measured  gravity  data.  The  gravity 
minimum  in  the  central  Alps  is  about  40  mGal  lower  in  our  map.  The  gravity  high  in  the 
Po  plain  is  well  approximated.  The  gravity  fit  is  worse  in  Austria  as  compared  to  detailed 
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Fig.  3.8a.  Gravity  effect  [mGal]  of  the  variation  of  the  seismic  basement  due  to  a 
density  contrast  of  -0.176  g/cm3,  in  relationship  to  a  normal  seismic 
basement  of  5.5  km  depth. 
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Fig  3.8b.  Gravity  effect  [mGal]  of  the  variation  of  the  Moho  due  to  a  density 
contrast  of  -0.278  g/cm3,  in  relationship  to  a  normal  crustal  thickness  of 
30.0  km. 
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Fig.  3.8c.  Gravity  effect  [mGal]  of  the  variation  of  the  lithosphere/asthenosphere 
boundary  due  to  a  density  contrast  of  +0.029  g/cm3,  in  relationship  to  a 
normal  lithospheric  thickness  of  105  km. 
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gravity  maps  (e.g.,  STEINHAUSER,  PUSTIZEK  1987)]  the  reason  for  this  is  probably  the 
poorer  crustal  model  in  this  region.  To  the  north  in  the  foreland  of  the  Alps  the 
interpretation  of  PLAUMANN  (1987)  is  confirmed. 

The  greatest  contribution  to  the  total  gravity  effect  is  that  of  the  Moho  variation.  It  is  the 
exclusive  cause  of  the  Alpine  gravity  minimum.  This  confirms  the  statistically  inferred 
aspect  of  the  good  isostatic  compensation  of  the  Alpine  structure.  The  mismatch  of 
40  mGal  in  the  axial  region  cannot  be  explained  by  changing  the  Mohorovicic  discontinuity 
within  the  limits  of  confidence  (Figs.  8.14,  8.15).  A  rock  layer  of  1  km  thickness  at  40  km 
depth  with  a  density  contrast  of  0.3  g/cm3  will  have  a  gravity  effect  of  less  than  10  mGal  at 
the  surface  and  a  change  of  Moho  depth  within  the  permissible  limits  will  not  explain  more 
than  50  %  of  the  residuum. 

Fitting  the  geoid  data  (see  Fig.  3.10)  in  the  same  manner  how  it  was  done  with  the  mean 
Bouguer  anomalies  leads  to  the  density  contrasts  represented  in  Table  3.5.  The  density  of 
the  topographic  masses  was  fixed  to  2.67  g/cm3  as  mentioned  above. 


Table  3.5.  Density  contrasts  at  the  layer  boundaries  from  three-dimensional  model 
computations  (fitting  of  geoid). 


Model 

Co 

c, 

c2 

c3 

R 

STDV 

4 

41.32  ±0.69 

-0.234  ±0.005 

m 

3.315 

5 

67.42  ±0.75 

-0.453  ±0.010 

-0.342  ±0.004 

0.88 

2.369 

6 

80.19  ±0.70 

-0.421  ±0.008 

-0.197  ±0.005 

-0.071  ±0.002 

1.898 

Co 

C,  ... 

c2  ... 
c3  ... 

R 

STDV  ... 


constant  of  fitting  [mGal], 

density  contrast  at  the  seismic  basement  [g/cm3], 
density  contrast  at  the  Moho  fg/cm3], 

density  contrast  at  the  lithospnere  /  asthenosphere  boundary 
[g/cm3], 

correlation  coefficient 

standard  deviation  of  the  model  [mGal], 


Model  4  ...  only  Moho  included, 

Model  5  ...  Moho  and  the  seismic  basement  included, 

Model  6  ...  Moho,  basement,  and  lithosphere/asthenosphere  boundary. 


The  density  contrasts  relative  to  the  sub-lithospheric  mantle  and  the  absolute  values  (again 
with  the  assumption  of  3.3  g/cm3  for  the  lithosphere)  are  as  follows:  lithosphere:  -  0.071 
(3.30)  g/cm3;  lower  crust  above  the  Moho:  -  0.268  (3.103)  g/cm3;  and  the  upper  crust 
above  the  seismic  basement:  -  0.689  (2.682)  g/cm3.  The  densities  derived  from  geoid  inver- 
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sion  are  also  quite  reliable  but  they  differ,  however,  from  those  derived  from  gravity  inver¬ 
sion  mainly  in  two  points.  First  the  significant  negative  density  contrast  at  the  lithosphere/ 
asthenosphere  boundary  has  to  be  noticed,  secondly  the  strong  negative  density  contrast  at 
the  seismic  basement.  In  Table  3.5  one  could  see,  that  the  densities  in  the  different  models 
are  even  more  dependent  on  the  number  of  influencing  boundaries. 

The  effects  of  topography  and  the  boundaries  with  the  density  contrasts  of  model  6  and  the 
total  geoid  effect  resulting  from  the  best  fit  of  the  input  data  are  shown  in  Figs.  3.11  a,b,c,d 
and  3.12.  The  greatest  contribution  to  the  calculated  geoid  is  that  of  the  topography;  the 
variation  of  the  Moho  has  not  the  strong  effect  on  geoid  as  on  gravity.  Although  the  geoid  is 
well  fitted  (comparison  of  Figs.  8.17  and  8.19  and  the  correlation  coefficient  of  0.92 
corroborate  this  fact)  the  density  contrasts  found  do  not  lead  to  a  linear  velocity-density 
systematics.  About  the  reasons  one  can  only  speculate.  Our  opinion  is  that  edge  effects  of 
the  seismic  model  are  not  entirely  eliminated  and  that  the  geoid  data  resolution  calculated 
from  a  spherical  harmonic  expansion  up  to  order  and  degree  360  is  too  coarse  to  include  the 
effect  of  such  small  areas  like  the  Alpine  lithospheric  roots  with  its  positive  density 
contrast.  Such  geoid  data  are  not  suitable  to  deduce  velocity-density  systematics,  but  its 
inversion  could  be  used  to  consider  the  results  of  gravity  inversion.  If  the  seismic  model  is 
good,  the  geoid  inversion  should  give  results  of  the  same  order  of  magnitude.  This  is  true  in 
case  of  our  model. 


3.6  Velocity-density  systematics 

The  term  "systematics"  for  the  relationship  between  density  and  seismic  velocity  is  used  to 
indicate  that  there  is  only  a  trend  with  much  real  scatter.  Solid  state  theory  cannot  yet 
establish  unequivocal  relationships.  Velocity  and  density  are  influenced  by  the  material  as 
such  (composition,  mean  atomic  weight)  and  by  the  packing  and  bonding  of  the  atoms 
(crystallography,  mineralogy),  i.e.,  by  many  parameters.  Thus  in  order  to  connect  gravity 
and  seismic  data  via  density-velocity  systematics  one  tries  to  establish  them  empirically. 
Mostly  one  refers  to  laboratory  results,  e.g.,  by  BIRCH  (I960,  1961). 

It  is,  however,  also  possible  to  take  the  reversed  route.  If  we  compare  the  density  contrasts 
found  with  the  two-  and  three-dimensional  gravity  models  (density  contrasts  found  from 
geoid  inversion  are  excluded)  and  the  average  seismic  velocities  vp  in  the  layers  of  our 
models,  we  find  a  nearly  linear  relationship  (Fig.  3.13).  The  density  contrasts  in  this  case 
have  been  transformed  to  absolute  densities  with  the  assumption  of  a  value  of  3.3  g/cm3 
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F’lg  3.10.  Geoid  undulation  [m]  calculated  from  a  spherical  harmonical  expansion  up 
to  order  and  degree  360  in  the  Alpine  area 
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Fig.  3.11b. 


Geoid  effect  [m]  of  the  variation  of  the  seismic  basement  due  to  a 
density  contrast  of  -0.421  g/cm3  in  relationship  to  a  normal  seis¬ 
mic  basement  of  5.5  km  depth. 
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Fig.  3.1  Id.  Geoid  effect  [m]  of  the  variation  of  the  lithosphere/asthenosphere 

boundary  due  to  a  density  contrast  of  -0.071  g/cm*  in  relationship 
to  a  normal  lithospheric  thickness  of  105  km. 
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within  the  uppermost  mantle.  We  find  the  coefficient  b  =  dvp/dp  to  be  about  5.4  ±  0.4 
[km/s]/[g/cm3].  This  result  is  relatively  high  as  compared  to  other  investigations.  For 
example,  JACOBY  (1975)  -  Fig.  3.15  -  found  values  about  4.8  ±  0.8  [km/s]/[g/cm3]  for 
long  range  models  or  profiles  in  North  and  South  America.  Laboratory  investigations 
(BIRCH  1960,  1961)  have  given  somewhat  lower  values  of,  for  example,  about  3.5 
[km/s]/[g/cm3]  for  basalt-type  rocks.  We  therefore  shall  not  overemphasize  our  linear 
relationship. 


Fig.  3.13.  Velocity-density  systematics  from  two-  and  three-dimensional  gravity 
inversion;  triangles:  two-dimensional  models;  asterixes:  three-dimensional 
models;  quadrangles:  three-dimensional  geoid  model  (for  comparison),  the 
gradient  dvp/dp  =  5.39  [km/s]/[g/cm3] 


Fig.  3.14.  Experimental  velocity-density  systematics  after  WOOLLARD  (1975). 
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However,  the  gravimetiic  results  as  ours  principally  suffer  from  their  insensitivity  to  purely 
depth-dependent  density  variations  in  the  earth.  It  is  quite  likely  that  there  is  generally 
such  a  component  of  purely  vertical  density  increase  in  the  real  density  variation  of  the 
crust  and  upper  mantle;  this  component  may  differ  from  region  to  region  and  thus  explains 
the  deviation  of  our  result  from  others  obtained  by  similar  methods  and  from  laboratory 
measurements. 
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Fig.  3.15.  Velocity-density  systematics  from  gravity  studies  after  JACOBY  (1915). 


This  is  no  real  handicap  for  the  present  investigation  and  for  the  aim  of  regional  gravity  pre¬ 
diction,  since  the  component  of  vertical  density  variation  does  not  influence  the  gravity  varia¬ 
tion  in  a  region  where  this  component  is  approximately  constant.  If  it  changes  on  larger  re¬ 
gional  scale  the  corresponding  gravity  variations  will  be  of  about  the  same  length  scale.  If  the 
studies  are  applied  to  large  regions,  the  vertical  density  variation  common  to  such  regions 
may  be  smaller.  If  smaller  and  larger  scale  studies  are  combined  it  should  be  possible  to  link 
the  gravity  prediction  scheme  to  the  scale  ( wavelength ,  harmonics )  of  gravity  field  representa¬ 
tions  of  sufficient  reliability.  In  any  case,  the  approach  to  the  question  of  gravity  prediction 
as  proposed  here  is  probably  best  for  scales  of  order  of  1 02  -  I03  km  and  not  good  on  very 
large  scales  (>  10 3  km).  If  the  seismic  data  are  good  and  the  gravimetric  data  are  good  m 
parf,s  or  for  relatively  long  wavelength  two-  and  three-dimensional  modelling  as  presented  in 
this  chapter,  we  believe  we  can  advance  to  a  finer  spatial  resolution  of  the  gravity  field  with  a 
reliability  comparable  to  that  of  seismic  data.  The  method  can  also  be  used  to  interpolate  or 
to  some  extent  even  to  extrapolate  the  gravity  field. 
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4.  A  BRIEF  SURVEY  OF  PROPOSED  (GEODETIC!  APPROACHES  IN  THE  PAST 


Special  emphasis  in  this  brief  survey  is  put  on  approaches  relevant  for  geodesy.  It  may  be 
by  no  means  complete.  The  idea  of  using  geophysical  information  when  predicting  gravity 
anomalies  is  not  new,  see  ( WILCOX  1974).  It  is  obvious,  that  the  solution  of  the  "inverse 
problem"  plays  an  important  role  in  geophysics.  A  mathematically-oriented  overview  of 
possible  methods  can  be  found,  e.g.,  in  TARANTOLA  (1987). 

Consequently,  we  may  categorize  the  work  done  up  till  now  into  three  groups: 

(1)  Geophysical  treatment  of  inverse  problem  theory.  As  mentioned  before,  one  can  find 
chapters  on  inverse  problem  theory,  e.g.,  the  determination  of  the  internal  structure 
of  the  crust  from  other  data  like  seismic  and  gravity  measurements  in  every 
textbook.  From  the  large  variety  we  only  like  to  mention  BOTT  (1971),  BULLEN 
(1975),  MEISSNER  (1986),  or  the  handbook  of  LANDOLT-BORNSTEIN  (1982). 
Recently,  the  complex  relationships  between  mass  heterogeneities  and  surface 
observables  such  as  topography,  gravity,  etc.  were  treated  by  various  authors  (see 
e.g.,  FLEITOUT,  FROIDEVAUX  1982;  KHAN  1977;  VIGNY,  FROIDEVAUX 
1988;  WOODHOUSE,  DZIEWONSKI  1984).  Especially  in  seismics,  the  term  "317- 
tomography”  is  now  used  with  respect  to  studies  of  the  earth’s  internal  structure. 
Quite  recently,  MORITZ  (1989)  presented  a  (geodetic)  approach  to  the  gravimetric 
inverse  problem. 

(2)  Indirect  use  of  geophysical  data  in  geodetic  methods.  Since  the  determination  of  the 
gravity  potential  in  geodesy  is  considered  to  be  a  boundary  value  problem,  where  the 
boundary  data  like  gravity  anomalies  etc.,  have  to  be  given  on  the  geoid,  it  is 
necessary  to  reduce  the  observations  from  the  physical  surface  of  the  earth  to  the 
geoid.  For  gravity  field  interpolation  a  certain  smoothing  is  anticipated  in  order  to 
reduce  the  prediction  errors. 

For  these  purposes  a  simple  density  model,  eventually  in  connection  with  isostatic 
compensation  hypotheses  and  a  digital  terrain  model  are  used  to  take  into  account 
some  kind  of  geophysical  evidence.  Thereby  the  density  between  the  surface  of  the 
earth  and  the  Mohorovicic  discontinuity  is  simplified  as  global  constant,  p0  =  2.67 
g  cm-3.  In  the  frequently-used  Airy-Heiskanen  isostasy  model  (see,  e.g.,  HEISKA- 
NEN,  MORITZ  1967,  p.  135f.)  the  mountains  of  constant  density  p0  float  on  a  dense 
underlayer  of  constant  density  p\  =  3.17  g  cm*3.  The  condition  of  floating  equilibrium 
is  (in  flat  earth  approximation) 
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t  A p  =  h  p0  (4-1) 

where  A p  =  px-  pQ,  h  is  the  topographic  height.  The  thickness  of  the  corresponding 
root,  starting  at  the  compensation  depth  T,  is  defined  by  t.  Assuming  standard 
values  (Ap  =  0.6  g/cm3;  p0  =  2.67  g/cm3)  leads  to  the  "rule  of  thumb"  t  =  4.45  h. 
The  crustal  thickness  under  mountains  is  T  +  h  +  t  ,  where  the  normal  thickness  (— * 
Moho  depth,  to  some  extend)  is  generally  assumed  to  be  T  =  30  km.  These  short 
descriptions  illustrate  more  or  less  the  (poor)  consideration  of  geophysical  data  in 
geodetic  methods. 

The  recent  developments  of  digital  (surface)  density  models,  see  for  example 
WALACH  (1987),  and  also  the  construction  of  Moho  models  (see  chapter  8.2.6)  will 
lead  to  a  significant  improvement  in  near  future  in  the  indirect  use  of  geophysical 
quantities  in  geodetic  algorithms  (SUNKEL  1986 ;  TSCHERNING  1979).  High- 
resolution  global  isostatic  earth  models  were  presented  by  SUNKEL  (1986),  see  also 
RJJMMEL  et  a 1  (1988).  More  detailed  investigations  on  "experimental  isostasy"  can 
be  found  in  ( DORMAN ,  LEWIS  1970 ;  LEWIS,  DORMAN  1970). 

In  a  more  general  sense,  one  can  understand  those  smoothing  approaches  with  respect 
to  mixed  models  (or  collocation)  also  in  such  a  way,  that  the  more  the  deterministic 
part  is  increasing  the  stochastic  treatment  of  the  rest  will  become  smaller  and 
smaller. 

Quite  recently,  GRAFAREND  (1989)  proposed  to  replace  *he  standard  topographic 
mass  reduction  by  a  volume  integral  over  the  datum-dependent  internal  Neumann 
kernel  and  anomalous  masses  6p  plus  a  surface  integral  over  the  surface-reduced 
internal  Neumann  kernel  and  the  anomalous  boundary  data  of  gravimetric  type. 

(3)  Direct  use  of  geophysical  data  in  geodetic  computations.  Direct  use  in  this  context 
means  that  we  (i)  either  use  the  geophysical  data  as  some  kind  of  (pseudo-) 
observations  in  our  geodetic  algorithms,  or  (ii)  set  it  into  correlation/regression  to 
gravity  field  functionals  like  geoidal  heights,  etc.,  or  (iii)  use  it  in  the  construction  of 
earth  gravity  models  consistent  with  internal  mass  distributions.  The  last  is  often 
also  a  prerequisite  for  methods  (i). 

With  regard  to  (i)  attempts  were  made  to  incorporate  density  into  a  stochastic 
concept  using  the  theory  of  random  processes  (JORDAN  1978,  MORITZ  1977, 
TSCHERNING  1976,  1977).  SANSO  (1980)  proposed  the  application  of  the  colloca¬ 
tion  principle  to  the  internal  densities  (internal  collocation).  Later  on  SANSO, 
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TSCHERNING  (1982)  represented  the  anomalous  gravity  field  by  a  combination  of 
an  anomalous  potential  Ts  harmonic  down  to  some  Bjerhammar  sphere  S0  and  an 
anomalous  potential  Tt)  generated  by  a  layer  of  topographic  masses  between  the 
sphere  S0  and  the  topographic  surface  S  (mixed  collocation).  Recently,  the  so-called 
"quasi-harmonic  inversion "  of  gravity  field  data  was  tested  ( TSCHERNING ,  STRY- 
KOWSKI  1987;  HEIN  et  al  1988).  The  central  problem  in  those  approaches,  the 
choice  of  norm  for  the  density  distribution  -  earlier  treated  by  KRARUP  (1970), 
SANSO  et  al  (1986),  TSCHERNING  (1974)  -  was  solved  in  such  a  way  that  density 
6p  is  supposed  to  be  a  solution  of  a  partial  differential  equation,  A(<5p/p(r))  where 
p(r)  is  a  polynomial  only  dependent  on  the  distance  from  the  origin  r,  and  A  is  the 
well-known  Laplace-0 perator.  In  addition,  indicator  functions  with  overlapping  sup¬ 
port  were  proposed  as  covariance  functions  ( HEIN  et  al  1988;  TSCHERNING  1989). 
VASSILIOU,  SCHWARZ  (1987)  discussed  different  methods  out  of  the  variety  of 
the  geophysical  literature  for  the  combination  of  gravity  with  other  geophysical  data 
for  the  solution  of  the  inverse  problem.  STRYKOWSKI  (1989)  presented  density 
autocovariance  functions  from  North  Sea  density  logs. 

With  respect  to  category  (ii)  example-wise  we  like  to  mention  studies  of  COLIC  et  al 
(1988),  FOTIOU  et  al  (1988),  GEISS  (1987),  LAMBECK  (1976),  McNUTT  (1980), 
RICARD  et  al  (1984). 

The  construction  of  reference  (normal)  density  models  consistent  with  the  outer  gravity 
potential  and  vice  versa  is  discussed  in  chap.  6.4  (MARTINEC,  PEC  1986a, b]  MATYSKA 
1987 ;  MESHCHERYAKOV  et  al  1986]  PEC,  MARTINEC  1984,  1988]  TSCHERNING, 
SUNKEL  1981). 
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5.  THE  INTEGRATED  GEODESY  ADJUSTMENT  MODEL 


As  mentioned  before,  our  goal  is  to  incorporate  the  data  types  density,  seismic  velocities 
and  displacements  in  the  direct  way  into  geodetic  algorithms  for  the  approximation  of  the 
gravity  Geld  outside  the  earth.  For  the  sake  of  an  unified  theory  we  will  try  to  derive 
(pseudo-)  observation  equations  fitting  into  the  general  scheme  of  integrated  geodesy 
adjustment.  The  mixed  model  approach  is  therefore  also  suited,  since  precise  relationship 
between  density  and  gravity,  for  example,  can  often  only  be  found  locally.  Necessary 
generalizations  for  regional  areas  imply  that  besides  such  an  (attempt  of  a)  deterministic 
mathematical/physical  relation  (more  or  less  local)  residuals  may  remain  which  can  be 
considered  as  forming  a  random  process.  This  is  in  correspondence  also  with  geophysical 
proposals  and  investigations,  to  consider,  for  example,  seismic  wave  propagation  and  densi¬ 
ty  perturbations  in  a  random  medium  ( CHERNOV  1960,  KORN  1987). 

We  review  here  the  general  principle  of  integrated  geodesy  ( EEG ,  KRARUP  1973,  HEIN 
1986).  Every  geodetic  and  geophysical  measurement  1  can  be  expressed  as  a  nonlinear 
functional  depending  on  one  or  several  position  vectors  x  =  (x,y,z)  in  space  and  on  the 
gravity  field  of  the  earth,  symbolically  written 

1  =  F(x,W)  (5-1) 

where  W  is  the  gravity  potential 

W  =  V  +  w2(x2  +  y 2)/2  .  (5-2) 

V  is  the  potential  of  the  gravitational  force  and  u)  is  the  angular  velocity  of  the  earth’s 
rotation.  By  (x,y,z)  or  (x1,X2,X3)  we  denote  the  Cartesian  coordinates  of  the  geocentric 
system  E  defined  as  follows:  the  origin  is  at  the  earth’s  center  of  mass,  the  z-axis  coincides 
with  the  (mean)  rotation  axis,  and  the  x-axis  goes  through  the  (mean)  Greenwich  meridian. 
We  presume  further  that  the  observations  are  corrected  for  time-dependent  geodynamic 
effects. 

The  scope  of  geodesy  is  now  to  determine  the  coordinates  of  material  points  on  the  surface 
of  the  earth  (and  in  space)  and  the  gravity  potential  including  its  functionals  by  the 
relation  (5-1).  Since  our  measurements  are  nonlinear  functionals  we  have  to  introduce 
approximate  values  for  the  linearization  process. 
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X  =  X0  +  6x 

W(x)  =  U(x)  +  T(x) 


(5-3) 

(6-4) 


xo  is  the  approximate  position  vector,  U(x)  is  some  kind  of  trial  value  for  the  actual  gravity 
potential  at  x  .  This  could  be  a  so-called  normal  potential  belonging  to  an  arbitrary 
ellipsoid  or  any  other  reference  potential,  e.g.,  a  low-order  harmonics  expansion  derived 
from  satellite  geodesy.  For  the  sake  of  comparison  of  the  results  with  those  of  other 
computations  one  can  use  one  of  the  adopted  reference  systems  of  the  International 
Association  of  Geodesy.  As  done  hitherto  we  will  call  T(x)  the  disturbing  potential  and  we 
will  assume  further  that  both,  T(x)  and  6x  =  (Sx,fy,Sz)r  =  (&C!,&C2,&C3)T  are  small 
quantities,  which  allow  us  to  work  with  a  linear  model. 

Thus  we  get  for  (5-1) 


I  =  F(x0+  6x,  U  +  T)  .  (5-5) 

Applying  Taylor’s  theorem  at  xo  =  [x?,x8,x§]T  and  restricting  ourselves  to  first-order  terms 
in  the  series  it  follows 

3 

1  =  F(x0,U)  +  £  Fx.  (x0,U)  6x{  +  L(T)  (5-6) 

i  *  1 


where  L(T)  is  a  linear  operator  applied  on  T  .  By  the  substitutions 

a  =  l-F(x0,U)  (5-7) 

ai  =  ^t(*o.U)  (5-8) 


we  get  the  general  linear  observation  equation  of  the  form 

=  aT  6x  +  L(T),  or  (5-9) 

considering  noise  n,  in  matrix  form  (A  and  R  are  corresponding  design  matrices), 

=  A<5x  +  Rt  +  n  (5-10) 
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where  the  first  expression  of  the  right  hand  side  represents  the  metric  coordinates  part  and 
the  second  the  (functionals  of  the)  disturbing  potential,  in  detail 


t  =  [T,  dT/dxu  PT/foi  ftcjjT  i j  =  {1,2,3}  . 


(5-11) 


Simple  mapping  and  scaling  of  the  quantities  in  (5-11)  yields 


-  geoidal  heights 


-  gravity  disturbances 


N  =  T/7,  (7  is  normal  gravity), 

*  _  STT 
^  "-^T> 


-  gravity  anomalies  Ag  =  -  ^  T 

(h  indicates  (approximately)  the  direction  of  the  normal  to  the  ellipsoid), 


-  deflections  of  the  vertical, 


north-south  component 


,  1  <9T 

^  K7  U<p 


east-west  component  r,  = 

where  R  is  the  mean  earth’s  radius,  and  {<p, A)  are  spherical  latitude  and  longitude, 


-  anomalous  gravity  gradients. 

The  observation  equation  system  (5-10)  can  be  solved  using  the  hybrid  minimum  norm 
condition 


nT  C„n  n  +  tT  Ktt  t  =  min 


(5—12) 


where  Cnn  and  Ku  are  the  covariance  matrices  of  n  and  t  .  The  covariance  matrix  Ktt  can 
be  derived  from  a  (global)  covariance  model. 


The  solution  of  a  general  collocation  model  of  type  (5-10)  can  be  found,  e.g.,  in  ( MORITZ 
1980,  p.  U If).  The  unknown  coordinates  are  given  by 


6x  =  (atd_1a)_1  AtD  1  a 


(5-13) 
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and  the  functionals  of  the  disturbing  potential  at  the  observation  points  are 


1  =  Ktt  Rt  -  A  6x)  .  (5-14) 

Since  in  the  collocation  model  an  interpolation  of  the  stochastic  part  is  implicitly  built  in, 
(5-14)  can  be  also  used  for  the  determination  of  any  other  functional  s  of  the  gravity 
disturbing  potential  at  any  other  station  when  knowing  the  corresponding  crosscovariances 
ILit 


S  =  KstETD  l(fl  -  Afic)  .  (5-15) 

The  error  statistics  are  given  by 

Exx  =  (ATD_1A)_1  (5-16) 

and 

Ess  =  Kss  -  KstD"'[I  -  A  (ArV~lA)-1  At  D'*]  Kts  (5-17) 

where  I  is  the  identity  matrix  and 

S  =  Cnn  +  R  Ktt  RT  ■  (5-18) 

For  the  estimation  process  above  there  are  certain  assumptions  necessary,  as 

E{n}  =  0  (5-19) 

E{t}  =  0  (5-20) 

E{s}  =  0  (5-21) 


where  E  =  £M  is  the  total  average.  E  is  the  expectation  operator,  and  M  describes  a 
homogeneous  and  isotropic  average  over  the  sphere. 

In  case,  that  we  assume  the  positions  x  to  be  known  where  our  measurements  are  taken,  the 
observation  equation  (5-10)  reduces  to 


fl  =  R  t  +  n  (5-22) 

and 

A  6x  =  0  (A  =  0)  in  (5-13)  to  (5-17)  . 
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Central  point  in  a  possible  inclusion  of  density  p  and  seismic  velocities  v  in  the  approach 
above  (consideration  in  t  eventually)  is  the  definition  of  the  necessary  covariance  matrix 
Ku  which  is  based  on  the  autocovariance  functions  cov(p,p)  and  cov(v,v)  as  well  as  the 
necessary  crosscovariances  cov(/?,t),  cov(v,t).  The  fulfillment  of  the  necessary  conditions 
(5-20),  (5-21)  assumes  that  density  as  well  as  seismic  velocities  are  considered  to  be 
random  which  they  are,  of  course,  in  reality  not!  But  for  the  application  of  the  theory  of 
random  processes  it  can  be  artificially  derived  through  trend  elimination;  in  other  words, 
subtracting  an  approximate  deterministic  model  leads  to  the  corresponding  anomalous 
(random)  quantities  6p,  6v  . 

For  the  deterministic  relations  underlying  the  crosscovariance  propagation  we  have  to  find 
some  (simple)  formulation  which,  obviously,  concerns  the  so-called  inverse  problem.  As  a 
consequence,  possible  constraints  have  to  be  considered  in  order  to  achieve  uniqueness  of 
the  solution.  However,  the  reader  should  be  reminded  to  the  discussion  in  chapter  2,  in 
particular,  see  Fig.  2.1. 

The  application  of  mixed  model  approaches  can  be  nowadays  found  in  many  scientific 
disciplines.  It  is  always  just  suited  when  the  deterministic  relationships  for  the  explanation 
of  mathematical/physical  phenomena  are  insufficiently  known,  so  that  no  other  choice  is 
than  to  live  with  large  model  discrepancies  or  to  treat  those  residuals  in  a  random  way. 
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6.  DENSITY  IN  AN  INTEGRATED  GEODETIC  APPROACH 


6.1  Minimum  norm  solutions 

In  the  following  we  like  to  discuss  (simple)  models  for  the  relationship  between  density  and 
the  gravity  potential  using  a  random  approach.  Some  of  them  are  already  applied  or  in  the 
context  of  this  work  numerically  tested. 


Generalized  Poisson  equation 

One  of  the  simplest  models  would  be  to  start  from  the  generalized  Poisson  equation 
( HEISKANEN ,  MORITZ  1967,  p.  47) 

AW  =  -4xGp  +  2wJ  (6-1) 

where  A  is  the  Laplace  operator,  W  the  gravity  potential,  G  the  gravitational  constant,  and 
u)  is  the  angular  velocity  of  the  earth’s  rotation.  The  potential  of  gravity  W  is  the  sum  of 
the  potential  of  the  gravitational  force  V  and  the  potential  of  the  centrifugal  force  4>. 


W  =  V  +  *  (6-2a) 

with 

$  =  0.5  d1  (x2  +  y2)  (6-2b) 

where  (x,y,z)  are  the  cartesian  coordinates  in  a  global  geocentric  system. 

We  linearize  by  introducing  (see  also  chap.  5) 

W  =  U  +  T  ,  (6-3) 

p  =  po  +  6p  .  (6-4) 

U  is  a  model  or  normal  potential  which  has  to  be  consistent  with  the  corresponding  model  or 
normal  density.  T  is  the  anomalous  (disturbing)  potential  and  bp  the  anomalous  density 
function.  We  choose  U  and  po  in  such  a  way  that  the  anomalous  quantities,  bp  and  T,  can 
be  considered  to  be  random  or  stochastic.  Inserting  (6-3),  (6-4)  into  (6-1)  and  solving  for 
bp  we  get 
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bp  =  —  (4  7r  G)'1  (AU  +  AT)  -  pa 


(6-5) 


where  the  anomalous  density  function  is  given  by 

3 

bp  =  —  (4  7r  G)'1  AT  =  -(4  ttG)'1^  32T/5xi  (6-6) 

i  =  1 

and  the  normal  density  is 

3 

Pa  =  —  (4  7r  G)'1  AU  =  -(4  tG)-1^  d2U/dx\  (6-7) 

i  =  1 

3  2  2 

Considering  (6-6)  as  a  linear  observation  equation  where  AT  =  £  d  T/dxj  is  part  of  the 

i  =  I 

signal  vector  t,  see  (5-11),  we  can  define  the  necessary  covariance  function  by 
crosscovariance  propagation  from  those  of  the  second-order  derivatives  of  the  potential,  for 
example, 

cov  (bp,  bp)  =  -  (4  7T  G)-2  cov  (AT,  AT)  (6-8) 

The  main  drawback,  however,  lies  in  the  fact  that  the  generalized  Poisson  equation  is  valid 
only  in  the  interior  of  the  earth,  but  we  are  interested  in  the  determination  of  the  gravity 
potential  outside  the  earth.  The  potential  inside  the  earth  is  not  harmonic.  Second-order 
derivatives  have  jumps  at  discontinuities  of  density.  Thus,  (6-6)  can  be  considered  only  as 
one  possible  solution  out  of  the  solution  space,  a  so-called  zero-potential  density  (MORITZ 
1989). 


Density  as  quasi-harmonic  function 

In  order  to  overcome  the  difficulties  mentioned  above  a  proposal  was  made 
(TSCHERNING,  STRYKOWSKI  1988)  to  consider  density  as  a  quasi-harmonic  function 
using  the  condition 

A  Wp(r))  =  6  (6-9) 
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p(r)  is  a  polynomial  only  dependent  on  the  radius  of  the  point  under  consideration. 
Assuming  p(r)  =  r«  and  considering  (6-4)  we  get  using  a  spherical  harmonic  expression, 

for  the  covariance  function  (or  reproducing  kernel)  of  the  disturbing  potential 

Ctt  (P,Q)  =  £^(^]i+lpi  (cos#)  (6-10) 

i  =  2 

with  the  degree  variances  o\ 


(6-11) 


C jj  ,  Sy  are  the  normalized  coefficients  of  the  spherical  harmonics  expansion,  Pi  are  the  Le¬ 
gendre’s  polynomials,  if  is  the  spherical  distance  between  points  P,  Q  .  GM  is  the  product 
of  the  gravitational  constant  and  the  mass  of  the  earth,  R  the  mean  radius  of  the  earth,  and 
r  =  r(P),  r'  =  r(Q)  are  the  radial  distances  of  P,  Q. 


For  the  covariance  function  of  the  (anomalous)  density  function  we  get 


^6p6p(P>Q)  _  £(-')"  Pi(cos^) 

i  *2 


(6-12) 


with 


2 

r  \ 


2 
O  i 


(rsr  tp) 


(2i  +  l)a 


(2i  + 
RS" 


3  +  n)a 


and  the  crosscovariance  function  between  the  disturbing  potential  T  and  6p  is 


(6-13) 


oo 

Ct6p(P,Q)  =  £.mpPi  (cost)  •  (6-14) 

i  =  2 

The  cross-  (or  auto-  )covariances  of  all  other  functionals  of  the  disturbing  potential,  Sj,  Sj  , 
can  be  derived  by  (the  inner  product) 
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cov  (Si,  Sj)  =  Si(P)  Sj(Q)  Ctt(P,Q)  • 


(6-15) 


As  easily  verified,  the  use  of  p(r)  =  rn  in  the  so-called  quasi-harmonic  inversion  results  in 
anomalous  density  values  which  attain  their  maxima  and  minima  at  the  boundary  of  the 
earth.  Density  values  at  deeper  localities  within  the  earth  are  just  a  smooth  mirror  of  the 
corresponding  situation  at  the  earth’s  surface.  This  is  in  contradiction  of  the  physical  reality 
where  density  increases  with  depth. 


Density  as  harmonic  function 

Quite  recently  MORITZ  (1989)  has  given  a  general  solution  for  the  (global)  gravitational 
inverse  problem.  A  general  set  of  continuous  density  distributions  for  the  sphere  consistent 
with  the  gravitational  potential  outside  the  earth  was  represented  using  radius-dependent 
polynomials,  spherical  harmonics,  and  generalized  matrix  inverses.  Based  on  geophysical 
evidence  appropriate  radius-dependent  functions  may  be  found  which  allow  to  determine 
the  anomalous  density  (after  subtracting  a  reference  density  model,  see  chap.  6.4)  within 
the  earth.  Although  still  much  more  work  has  to  be  done  in  applying  this  model  in  numeri¬ 
cal  studies,  it  seems  that  the  approach  is  suited  for  the  construction  of  a  global  density- 
gravity  covariance  model  which  could  be  an  extension  of  the  well-known  Tscherning/Rapp 
spherical  harmonics  covariance  model  (TSCHBRNING,  RAPP  1974). 


Other  choice-  of  norm  proposals 

In  ( SANSO  et  al  1986)  different  minimum  norms  are  discussed  for  the  density  distribution 
of  the  earth.  The  generally  used  L2-norm  (5-12)  implies  further  that  blocks  of  constant 
density  are  uncorrelated.  If  base  functions  with  overlapping  support  are  used  this  problem 
could  be  overcome.  It  leads  to  a  reproducing  kernel 

C6p6p(P,Q)  =  Ii(P)  Ii(Q)  /  Vi  (6-16) 

where  I;  is  the  indicator  function  of  the  i-th  block  with  volume  v*  (TSCHERNING  1989). 
However,  the  drawback  there  is  that  the  density  covariance  function  has  values  equal  to  one 
and  zero  at  the  earth’s  surface  (HEIN  et  al  1988).  In  addition,  crosscovariance  propagation 
between  density  and  gravity  field  functionals  was  only  done  numerically  using  Newton’s 
attraction  integral  instead  of  (needed)  analytical  solutions  for  practical  applications. 
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Thus,  stronger  norms  or  combination  of  different  criteria  have  to  be  tried.  The  H1,J-norm  as 
proposed  by  SANSO  et  al  (1986)  seems  to  be  more  realistic 

Jra  6p2  +  b  | V(£p)| 2  dft  =  min  (6-17) 

ft 

since  it  can  take  into  account  lateral  variations  of  density  within  the  earth’s  crust  using  the 
horizontal  gradient  operator  V  in  (6-17). 

This  principle  can  be  extended  using  the  so-called  mixed  collocation  proposed  by  SANSO, 
TSCHERNING  (1982).  They  propose  to  split-up  the  gravity  disturbing  potential  into  a 
part  associated  with  an  internal  sphere  fto  of  the  earth  and  a  topographic  layer  including, 
for  example,  all  masses  between  the  topographic  surface  and  the  Moho.  At  the  first  part  the 
L2-norm  could  be  applied  whereas  to  the  second  a  Sobolev-type  of  norm  like  in  (6-17)  could 
be  tried, 

f  a  6p2  dft0  +  fb  |V(^)|2d(ft-ft0)  =  min.  (6-18) 

ft0 


6.2  The  attenuated  white  noise  statistical  gravity  model 

A  self  consistent  covariance  model  was  developed  1976  by  Heller  (HELLER,  JORDAN 
1979;  JORDAN  1978;  JORDAN,  HELLER  1978)  which  is  able  to  consider  also  topography 
and  density  contrasts  within  the  earth. 

The  spherical  model  for  the  autocovariance  function  of  the  disturbing  potential  T  has  the 
form 


CTT(r,  r.i,)  =  P’.Ct-P/jPC.T  _  (6-19) 

[  R*— (R-D/2 )  *1  [r2r'2+(R-D/2)*-2(R-D/2)2  r  r'  cos li'J 15 

where  ip  =  |  P— Q  |  is  the  spherical  distance  between  two  points  P,Q  on  the  earth’s  surface, 
R  is  the  mean  radius  of  the  earth,  and  r  =  r(P),  r'  =  r(Q)  are  the  radii  of  P,Q.  The  two 
free  parameters  (besides  ip)  scaling  the  model  are  the  variance  C0t  of  the  disturbing 
potential  and  the  characteristic  depth  D. 
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Model  (6-19)  was  derived  from  Poisson’s  integral  for  upward  continuation,  assuming  that 
the  disturbing  potential  is  uncorrelated  for  any  arbitrary  small  neighbourhood  on  the 
surface,  so  that  its  white  noise  process  distribution  looks  like 


Ctt  {t) 


A  S(t) 
2ttR2  s  i  nip' 


(6-20) 


where  6(ip')  is  the  Dirac  delta  function  and  A  is  the  spectral  density  of  the  impulse. 
Because  of  this  property  (6-20)  the  name  ’’attenuated  white  noise"  was  given  to  the  model 
(6-19),  since  the  surface  integral  of  (6-20)  over  the  sphere  with  radius  R  yields  the  value  A. 


For  local  applications  plane  approximations  of  (6-19)  are  useful.  Thus,  the  corresponding 
asymptotic  relation  is  obtained  by  considering  R  -»  od  or  (D/R)  -*  0,  ip-*  0  in  (6-19)  resulting 
in 


C  tt(u,v,z,z'  ) 


4D2  (  2D+z  +  z  Q  CoT 
[u2+v2+(2D+z+z')2j  15 


(6-21) 


where 

u  =  x'-x  =  R^cosa  (6-22) 

v  =  y'  -y  -  Resina  (6-23) 

The  points  under  consideration  have  the  three-dimensional  orthogonal  plane  coordinates 
P(x,y,z  Q(x',y',z')  which  are  related  by  (6-22),  (6-23)  to  their  spherical  counterpart,  a 
is  the  azimuth. 

JORDAN  (1978)  has  also  proposed  to  superimpose  many  models  of  type  (6-19)  for  more  or 
less  global  applications  each  representing  an  additional  white  noise  shell  (layer)  of  the 
earth’s  interior.  This  explains  also  why  he  denotes  bp  as  density  contrasts.  As  long  as  we  are 
working  only  with  one  shell  (layer)  the  density  contrast  between  outside  the  earth  and  the 
Grst  (surface)  density  layer  corresponds  to  anomalous  density  bp  used  throughout  this  re¬ 
port. 

The  auto-  and  crosscovariances  of  the  other  gravity  field  functionals  can  be  derived  by 
common  covariance  propagation  using  the  well-known  potential  relations  between  them. 
For  the  reader  who  is  interested  in  them,  we  refer  to  Appendix  A. 
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Interesting  in  this  context  is  how  the  link  to  density  (contrasts)  is  established.  It  is  founded 
by  the  theorem  of  Chasles  (HEISKANEN,  MORITZ  1967,  p.  IS)  which  states  that  the 
potential  V,  at  any  point  P  outside  the  equipotential  surface  S,  of  a  surface  layer  of  density 

p  =  -  (4xG)  >  dV/dn  ,  (6-24) 

is  the  same  as  that  of  the  attracting  body  itself,  n  is  the  outer  normal.  Thus,  the 
corresponding  relation  with  respect  to  the  density  contrast  bp  and  the  disturbing  potential 
T,  at  a  point  P(0,A,r),  with  spherical  coordinates  co-latitude  0,  longitude  A,  radius  r, 

Sp(0, A,r)  =  -  (4*G)-‘  T(0,A,r)  6( r  -  R*)  (6-25) 

where 

R*  =  R-r.  (6-26) 

By  6  the  Dirac  delta  function  is  denoted,  and  r  is  the  depth  of  the  considered  density 
contrasts  (in  the  case  of  many  layers).  The  derivative  in  (6-25)  is  evaluated  at  Rt  =  R*  . 

The  density  contrast  autocovariance  function  corresponding  to  (6-25)  and  the  autocovari¬ 
ance  of  the  disturbing  potential  Ctt  is 

<Wr,r',tf<)  =  (47tG)~2  CTT  •  #r  -  R*)  6{ r'  -  R*)  (6-27) 

For  gravity  field  determination  using  density  information  the  crosscovariance  between  free- 
air  gravity  anomalies  AgF  and  density  contrasts  6p  necessary  to  apply  least-squares  predic¬ 
tion  is  definded  by 


-  (4'G)-'  $F  [  ir  +  f]  C'T  -  R*)  ■  (6-28) 

Note :  When  considering  in  the  expressions  above  only  one  density  layer,  the  inversion  of 
matrix  D  in  the  prediction  formula  (5-15)  becomes  instable,  in  the  asymptotic  form  even 
singular  in  case  that  gravity  anomalies  and  density  data  as  observations  are  introduced. 
From  appendix  A  we  can  easily  deduce  that  the  covariances  Cgp6p  and  CAggp  as  well  as 
CgpAg  anc*  ^AgAg  are  *'near  dependent  then: 

C6,)&P(Z.Z’.V’)  =  (4xG)'<  6{z)  CAg6p(z,z',^) 
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C6pAg(z’z’’^  =  (4xG)'‘  6(z)  CAgAg(z’z’’^ 

Thus,  the  processing  of  a  combined  data  set  is  not  possible.  One  possibility  to  overcome 
these  difficulties  when  considering  only  one  layer  could  be  the  replacement  of  the  Dirac 
delta  function  S(z)  by  a  function  decreasing  linearly  or  even  exponentially  with  increasing 
height  or  height  difference  (z’-z). 


As  mentioned  by  JORDAN  (1978)  the  model  above  is  appropriate  for  density  contrasts 
which  are  not  compensated  isostatically.  He  has  therefore  developed  also  regional  isostatic 
compensation  models  for  terrain  and  crust-mantle  contrasts  based  on  a  floating  elastic  crust 
and  its  mean  (numerical)  material  constants  (Lame’s  parameters,  etc.).  Starting  point  of 
this  derivation  is  an  equation  for  the  displacement  zm  of  the  Moho  under  loading  which 
looks  like  the  one  treated  by  TURCOTTE,  SCHUBERT  (1982,  p.  122).  The  resulting 
autocovariance  model  for  gravity  anomalies  considering  regional  isostatic  compensation  is 


CAgAg(r’A’^=  C8p  4x202 


(1—1.5  r2/A2)  2A4  (1—1.5  r  2/Ai) 

(l+r2/A2)  3  5  Af  (l+r2/A|)  3  5 


|  A4  (1— 1.5r  2/A]) 
A$  (l+r2/A§)  3  5 


(6-29) 


where  r  =  R  ip  .  A  is  a  free  parameter  (characteristic  distance)  in  (6-29)  which  has  to  be 
determined  from  the  empirical  covariance  function.  C6p  is  the  variance  of  the  terrain 
density  contrasts.  Further,  we  have 


A,  =  A  +  0.693  d  +  H  (6-29a) 

A  3  =  A  +  1.386  d  +  2H  (6-29b) 


where  H  is  the  average  thickness  of  the  crust,  and  d  is  an  elastic  parameter  which  can  be 
computed  from 


d4 


c 

g  ( Pm  ~  Pc. ) 


(6-29c) 


with 


c  =  (6'29d) 

where  A,  p.  denote  the  Lame’s  constants  of  the  crust,  pm  is  the  mantle  density,  and  p o  -  1.03 
[g  cm  3]  for  oceans,  po  =  0  for  continents. 
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6.3  A  new  approach  to  link  gravity  and  density 

As  the  reader  might  have  realized  from  the  preceding  chapters,  in  particular  6.1,  an  optimal 
approach  can  only  be  found  when  considering  also  the  physical  relationship  between  gravity 
and  density.  The  following  proposal  is  mainly  based  on  the  theory  of  the  earth’s  isostatic 
response  to  a  concentrated  load  as  used  in  ( DORMAN ,  LEWIS  1970;  LEWIS,  DORMAN 
1970). 

For  the  derivations  afterwards  we  recall,  that  the  Fourier  transform  of  f(x)  is  defined  by 

QD 

§[f]  =  F(s)  =  J  f(x)e2l'1<— '  dx  (6-30) 

-tD 

and  its  inverse  transform  by 

CD 

3-i[F]  =  f(x)  =  f  F(s)  eJl'1<— ’  ds  (6-31) 

“ QD 

where  x,  s  can  be  also  two-  or  three-dimensional  vectors,  and  consequently,  the  integration 
two-  or  threefold. 

Convolution  (denoted  by  *)  of  two  functions  f(x),  g(x)  is  defined  by 

(D 

h(x)  =  f(x)  *  g(x)  =  J  g(u)  •  f(x  -  u)  du  (6-32) 

-QD 

The  convolution  itself  is  also  a  function  of  x  . 

6.3.1  Newton’s  attraction  integral 

We  start  with  the  gravity-density  relationship  by  Newton’s  attraction  integral, 

W(x)  =  G  J p{ y)  l'1  dvE(y)  +  $(x)  (6-33) 

ve 
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where  W  is  the  gravity  potential,  G  is  the  gravitational  constant,  p  is  the  density,  vE  is  the 
volume  of  the  earth,  x  =  (xt,  x2,  X3)  and  y  =  (yi,  y2,  ys)  are  threedimensional  vectors,  and 
1  =  |  x  -  y[  is  their  distance.  The  potential  of  the  centrifugal  force  is  given  by 

4>(x)  =  0.5  u)2  r2  cosfy  ,  (6-34) 

where  u>  is  the  angular  velocity  of  the  earth’s  rotation,  r  =  r(x)  is  the  radius  of  a  spherical 
earth,  and  <p  =  <^(x),  A  =  A(x)  are  the  spherical  latitude  and  longitude,  resp.  As  already 
mentioned  we  decompose  density  (6—4)  in  a  reference  or  model  part  po  and  an  anomalous 
(irregular)  term  6p,  the  last  forming  a  scalar  stationary  random  process  with  zero  mean, 

p(x)  =  P o(x)  +  6p(x)  .  (6-35) 


For  example,  the  model  density  po  can  be  defined  by  a  laterally  homogeneous  function  as 
done  in  some  earth  density  models  (see  chap.  6.4), 

p(<p,  A,  r)  =  p0(r)  +  Sp{ip,  A,  r)  .  (6-36) 

This  leads  to  a  split-up  of  (6-33)  into 

W(x)  =  G  J po(y)  l'1  dvg(y)  +  G  f  6p(y)  H  dvE(y)  (6-37) 

vg  vg 

where  the  first  integral  of  the  ^ht  hand  side  has  to  be  consistent  with  the  normal  or  refer¬ 
ence  potential  U  (see  also  6.*;  including  4>,  and  the  second  one  with  the  disturbing  or 
anomalous  potential  T.  We  will  further  assume  that  the  kernel  of  the  second  part  behaves 
in  such  a  way  that  the  integration  has  to  be  carried  out  only  in  a  local  region  with  limited 
volume  v,  so  that 


T(x)  =  G  f  6p( y)l-'dv(y). 
v 


(6-38) 


6.3.2  The  isostatic  response  function  m  space-domain  representation 

The  gravity  disturbance  (or  also  gravity  anomaly  in  flat-earth  approximation)  is  the 
vertical  derivative  of  T 
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(6-39) 


fe(x)  =  fp(x)  =  G  f6p(y)6(x, y)dv(y) 

v 

where  <?(x,y)  is  a  Green’s  function 

«M)  =  IfO')  =  3f(ls-*|-')  (6-*0) 

and  r  is  the  earth’s  radius. 

In  the  further  approach  we  assume  (using  spherical  coordinates  tp,  A,  r)  that  the  change  in 
density  6p  (v>,A,r)  at  a  given  depth  under  a  point  due  to  isostatic  response  is  a  function  of 
depth  alone  and  location  alone, 

6p(tp,\, r)  =  JJ  h  (tp' ,\')  6p  (t,  tp  -  tp' ,  \  cosp  dtp' d\' 

Up  (tp, A,r)  =  h  *  6p  (r,  tp-  tp' ,  A  -  A')  —  h  *  6p  (r,  it)  (6-41) 

where  *  denotes  the  convolution  (6-32)  of  the  functions  based  on  the  coordinates  tp,  A.  ip  is 
here  the  central  angle  between  (tp,  A)  and  (tp' ,  A'),  defined  by 

cosV>  =  cosOco&O'  +  sin0sin0'cos  (A  -  A')  (6—42) 

with  0  =  90°  -  tp  Thus,  the  (isostatic)  density  changes  are  now  expressed  as  a  product  of  a 
function  of  depth  alone  and  function  of  position  alone.  6p(i,ip)  is  the  characteristic  density 
change  due  to  a  unit  load,  and  h  =  h(v>,A)  is  the  elevation. 

Considering  (6—41)  in  (6-39)  we  get  (in  spherical  approximation) 

tfg(vJ,A,r)  =  h  *  q(r,tp)  (6—43) 

where  q(r,V»)  is  the  isostatic  response  function, 


te 

q(r,^)  =  G  J  6p  (r',ip)  ♦  Q  (f  dr'  .  (6-44) 

o 

booking  more  in  detail  to  (6-43)  <5g  must  be  the  Bouguer  anomaly  Ago  minus  the  gravity 
effect  due  to  variations  in  density,  Agp, 

h  =  Ago  -  Agp  (6-45) 
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6.3.S  The  isostatic  response  function  in  frequency-domain  representation 


As  done  in  other  methods  recently  -  see,  e.g.,  terrain  correction  determination  -  we  like  to 
use  the  advantages  of  transforming  the  formulas  in  the  frequency  domain  by  means  of  the 
Fourier  transform.  We  will  assume  that  due  to  the  local  applications  having  in  mind  a 
plane  approximation  is  sufficient  (possibly  supplied  by  small  corrections). 

Considering  (6-45)  in  (6-43)  the  isostatic  response  function  in  frequency-domain 
representation  is  given  by  (6—43)  using  (6-30), 


Q  =  S[q] 


S[AgB]  -S[Agp] 

*IKJ 


(6-46) 


Since  Agp  is  the  gravity  effect  due  to  variations  of  density  it  can  be  assumed  that  the 
crosscorrelation  between  Agp  and  topography  tends  to  zero.  Thus,  the  real  part  of  the 
isostatic  response  function  in  frequency-domain  representation  is 


Q(IQI)  i  Re 


(6-17) 


6.3.4  The  power  spectral  model  for  gravity  anomalies  ( disturbances ) 

The  autocovariance  function  for  6g  forming  a  scalar  random  field  is  given  by 

C6r6r(x)  =  M  (6g(y)  6g(x  +  y)}  (6-A8) 

where  M  is  a  suitable  averaging  operator.  Analogously,  we  have  for  the  density  anomalies 

6p 

C6p«pOO  =  M  Wx)  6p(x  +  y)}  .  (6-49) 

Using  the  Wiener-Khinchine  theorem  we  can  express  the  covariance  functions  above  in 
terms  of  their  power  spectral  densities  Cgg6g,  Sgpgp 

ChRgR(x)  =  /  SgRfeK(s)e?,,,(-)  ds  (6-50) 

-w 
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(6-61) 


<W*>  = 

-OD 

where  Sgg8g,  S8p8p  are  defined  by  the  Fourier  transform  (6-30).  Of  course,  the  inverse 
relations,  expressing  the  power  spectra  in  terms  of  their  covariance  functions  hold  too,  using 
the  (inverse)  Fourier  transform  (6-31). 

We  are  starting  our  derivation  now  using  (6-48)  and  inserting  Newton’s  attraction  integral 
(6-33).  After  reordering  the  integration  we  have 

C6g6g(x)  =  G2f  f  dx  dy'  </(x,y)  0(x',X')  M  (6p(y)  6p(x  +  y')}  (6-52) 

V  V 

where  M {6p(y)  6p(x+y')}  =  C6p6p(x)  represents  the  autocovariance  function  of  the 

density  anomalies  (6-49).  We  now  insert  in  (6-52)  the  corresponding  relation  in  terms  of  its 


power  spectrum  (6-51)  and  get 

C6g6g(x)  =  G*ff  dx  dy'  £(x,y)  S(x',r)  /S8p8p(s)  e*"-'1'  ds  .  (6-53) 

R.3 

The  two-dimensional  power  spectral  density  of  (6-53) 

S8g8g(s)  =  9[C6g6g(x)]  (6-54) 

il  of  the  form  (*  means  here  conjugate  complex) 

S8g8g(80)  =  GJ  J S8p8p(s)  I(s',x,y)  I*  (s',x,y')  ds  (6-55) 

K3 

where 

I(s',x,y)  =  /5(x,y)  e  2l,‘!  ?dy  (6-56) 

V 

I(i',x,y')  =  feu,*')*2"-'1'  dy'  (6-57) 

V 
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and 


s'  —  (s'j,  s 2 ,  s3)  • 


Solving  the  integrals  (6-56),  (6-57)  analytically,  e.g.,  with  software  like  REDUCE,  and 
considering  that  the  special  Green’s  function  (6-40)  used  here  is  defined  by  (£(x,y') 
similarly) 


£(x,y)  =  - ^3  X3 ) - 

[(yi-xO2  +  (y7-x2)2  +  (y3-x3)2)0'5 


we  get  for  (6-56),  (6-57)  the  expressions 


(6-58) 


I 


( i-l)s  3 

S  1 2  +  +  S32 


(6-59) 


I4 


s?  +  s^2  +  s32- 


(6-60) 


Thus,  expression  (6-55)  represents  a  gravity  disturbance  power  spectral  model  as  function 
of  the  power  spectral  density  anomaly. 


6.S.5  The  gravity  disturbance  covariance  model 

As  always  done  in  geodesy  we  choose  an  analytical  model  for  an  autocovariance  function 
whose  free  parameters  are  determined  by  a  best  approximation  to  the  empirical  derived 
covariance  values.  In  this  case  we  take  the  well-known  Gauss’  function  to  describe  the 
auto  covariance  function  of  density  anomalies, 

C6„tp(x)  =  C8„(2»)-o.*e-ls’l/2a  (6-61) 

The  free  parameters  are  the  density  variance  Cgp  and  the  length  scale  parameter  a.  (6-61) 
has  the  power  spectral  density 

S8p6p(s)  =  C^-o.s  f  e- 1  — 2I /2a  e  -  dx  (6-62) 
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After  solving  the  integral  in  (6-62)  we  get 


SSpSp(a)  =  C|pa>1re-2’lalM+si+si)  (6-63) 

which  is  inserted  in  the  gravity  disturbance  power  spectral  model  (6-55)  using  (6-59), 
(6-60) 


SSg6s00  =  2GJ  C?p  t  /e 

R3 


-2x2a2s2 


spT 


Tp“T 


irJd- 


resulting  in  the  analytical  expression 


(6-64) 


S6g6g(8)  =  G2  C6°p  a  x 


(1  +  45r2a2(s/,2+S22)]  e  2,r2a2(si2+s22) 


16  (s  j2+Sj2)® 


0-5 


(6-65) 


Considering  the  Wiener-Khinchine  theorem  (6-50)  we  can  derive  from  (6-65)  the  covariance 
function  of  gravity  anomalies  (disturbances)  as  function  of  the  variances  Cgp  of  the  density 
autocovariance  function  C6p6p,  a,  and  the  distance.  Thus,  we  have  a  consistent  model  to 
predict  gravity  anomalies  (disturbances)  using  density  information  on  the  basis  of  the  theory 
of  stochastic  processes. 


6.3.6  Inferences  for  gravity  prediction  applications 

In  the  application  of  the  theory  mentioned  above  the  assumption  has  to  be  made  that  the 
regional  mechanism  of  compensation  is  of  the  same  or  similar  nature  in  the  whole  area 
under  consideration.  Several  approaches  are  now  possible  to  use  the  density  information  for 
gravity  field  approximation. 

Let  us  assume  that  after  subtracting  a  certain  (e  g.,  laterally  homogeneous)  density  refer¬ 
ence  model  anomalous  density  values  6p  are  available  forming  a  scalar  stationary  random 
process  with  zero  mean.  Then  we  are  able  to  dv  mine  the  empirical  anomalous  density 
autocovariance  function  (6-49)  which  can  then  be  fitted  to  the  corresponding  model  (6-61) 
using  adjustment  techniques  resulting  in  the  determined  free  parameters,  the  variance  CJp 
and  the  length  scale  factor  a.  Using  these  quantities  the  gravity  anomaly  power  spectrum 
(In  flat-earth  approximation),  see  (6-65),  is  determined.  By  applying  the  Wiener-Khinchine 
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theorem  (6-50)  the  gravity  autocovariance  function  is  defined.  In  a  similar  way  as  outlined 
in  6.34  we  can  also  derive  the  crosscovariance  cov{6g,6p).  Thus,  we  are  able  to  define  the 
necessary  covariance  matrix  Ctt  needed  to  apply  the  prediction  approach  described  in 
chapter  5  (possibly  in  combination  with  other  geodetic  /  geometric  observations). 

It  would  be  also  possible  to  use  the  isostatic  response  function  itself  as  gravity  predictor. 
The  spectral  form  (6^6)  could  be  either  determined  from  a  Bouguer  gravity  anomaly  field 
and  applied  in  unsurveyed  areas  showing  the  same  compensation  mechanism.  In  this  case 
the  Wiener-Khinchine  theorem  has  to  be  used  and  the  convolution  (6-43)  carried  out  with 
the  actual  topography.  If  density  information  is  available  in  the  area,  the  isostatic  response 
function  can  be  determined  numerically  by  (6-44),  in  addition,  resulting  most  likely  in  a 
more  improved  and  detailed  resolution  which  reflects  the  local  compensation. 

If  namely  the  compensation  is  assumed  to  be  due  to  a  change  of  density  as  a  function  of 
depth  Ap(r)  occuring  directly  beneath  the  load,  then  the  transfer  function  Q  (6-46)  can  be 
directly  constructed  by 


fE 

Q(s)  =  G  J  Ap(r)  e  2ltr-  dr  .  (6-66) 

o 

This  follows  from  applying  Hankel’s  transform  to  (6-44). 

Some  available  gravity  data  in  the  area  would  be  useful  to  transform  the  gravity  anomalies 
derived  via  (6-43)  into  the  right  frame  (adding  a  trend,  for  example). 


Estimating  the  isostatic  response  function  is  widely  applied  now  in  modern  geophysics.  A 
survey  is  found,  e  g.,  in  (LAMBECK  1988;  p-425).  The  corresponding  transfer  function 
(6-46)  can  be  used  to  estimate  the  mean  density  of  the  topography,  similar  to  the  Nettleton 
approach,  see  LEWIS,  DORMAN  (1970,  p.S878). 


(6-67) 

(6-68) 

(6-69) 
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where  Agp  are  free-air  gravity  anomalies,  Agt  are  the  terrain  corrections,  and  H  is  the 
elevation.  The  brackets  <  •  >  indicate  averaging.  The  transfer  function  Q'(s)  is  related  to 
the  gravity  field  of  the  compensation  and  topography. 

One  remark  to  the  type  of  available  density  data:  In  general  those  data  come  from  borehole 
investigations  and  characterize  only  the  parts  of  the  earth’s  crust  near  the  earth’s  surface. 
In  constructing  the  isostatic  response  function  according  to  (6-44)  the  integration  has  to  be 
carried  out  between  the  center  of  mass  and  the  surface.  Thus,  for  the  interior  of  the  earth 
certain  density  earth  models  have  to  be  used.  If  the  stochastic  approach  (see  chapter  5)  is 
applied  as  mentioned  before  those  models  can  be  used  also  to  achieve  the  zero  mean  of 
density  anomalies  as  required  for  the  application  of  the  theory  of  random  processes. 

On  the  other  hand  it  might  be  worthwhile  to  consider  the  possibility  of  constructing  always 
first  the  isostatic  response  function  (6-44)  in  order  to  use  it  as  some  kind  of  (pseudo- ) 
observation  in  least-squares  prediction.  The  advantage  of  using  it  instead  of  (anomalous) 
density  values  lies  in  the  fact  that  the  isostatic  response  is  a  function  depending  only  on  the 
horizontal  distance  ,  A-A'),  and  no  more  on  the  depth  (or  radius,  resp.)  as  it  is  the 

case  with  density.  Insofar,  such  a  quantity  fits  better  to  an  integrated  adjustment  of 
geodetic/geophysical  observations  using  least-squares  collocation  (or  prediction). 


6.4  Reference  (normal)  density  models 

In  chapter  6.1  -  see  (6-3)  or  (6-37),  resp.  -  the  reference  or  normal  potential  was 
introduced 

U(x)  =  G  /po(x)HdvE(x)  (6-70) 

ve 

where  po  is  the  so-called  normal  or  model  density.  Due  to  the  integration  of  different  type  of 
observations  in  the  least-squares  collocation  model  of  integrated  geodesy,  p0  has  to  be 
consistent  with  U,  e  g.,  U  =  f(po)-  Thus,  assuming  that  we  usually  work  with  a  normal  po¬ 
tential  model  adapted  by  the  International  Association  of  Geodesy,  the  corresponding  earth 
density  model  has,  in  principle,  to  come  from  the  inverse  gravimetric  problem,  p0  =  f(U), 
which  has  no  unique  solution  (MATYSKA  1987;  MORITZ  1989). 
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Thus,  the  reference  density  model  has  to  fulfill  the  following  conditions: 

(1)  The  mass  of  the  earth  M  has  to  be  reproduced  by  po 

M*J  Po(y)  dvE(y)  (6-71) 

ve 

(2)  The  normal  density  distribution  should  have  ellipsoidal  or  spherical  geometry 
consistent  with  the  normal  potential.  Analogously,  the  geometry  of  the  different 
layers  of  shells  should  be  of  the  same  form. 

Thus  we  have  (with  respect  to  an  ellipsoidal  normal  potential  of  the  earth)  to  fulfill 
the  relation 

=  C  ~M(ttB)/2  (6-72) 

where  J20  is  the  underlying  coefficient  of  the  spherical  harmonic  expansion  (dynamic 
flattening).  A,  B,  C  are  the  moments  of  inertia  defined  in  (HEISKANEN,  MORITZ 
1967,  p.  62) 


A  =  /  (y5  +  yl)  Po(x)  dv(y) 

(6-73) 

VE 

B  =  /(yi  +  y!)  po(y)  dv(y) 

(6-74) 

ve 

c  =  /(yJ  +  y22)  p0(x)  dv(y) 

(6-75) 

ve 

where  y  =  (yi,y2>y3)  are  geocentric  cartesian  coordinates. 

(3)  If  the  density  distribution  is  within  a  sphere  with  mean  radius  R 

R  =  (6-76) 
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and  the  reference  potential  generating  body  is  rotating  with  angular  velocity  u,  our 
density  earth  model  is  consistent  with  the  normal  potential,  a,  b  in  (6-76)  are  the 
two  axes  of  the  ellipsoid.  For  a  spherical  earth  model  we  have  to  set  a  =  b. 

It  might  be  possible,  in  particular  with  respect  to  very  local  applications,  that  a  simple 
geometric  trend  function  can  be  used,  too,  in  order  to  generate  our  anomalous  densities  6p. 

Various  models  and  proposals  are  already  developed  in  the  geodetic  /  geophysical  literature. 
The  reader  is  referred  to  BULLEN  (1975),  KHAN  (1982),  MARTINEC,  PEC  (1986  a,b), 
MESHCHERYAKOV  et  al  (1986),  MORITZ  (1989),  PEG,  MARTINEC  (1984),  RUMMEL 
et  a 1  (1988),  SVNKEL  (1986),  TSCHERNING,  SUNKEL  (1981). 
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7.  SEISMIC  AND  SEISMOLOGICAL  DATA  IN  AN  INTEGRATED  APPROACH 


In  the  following  a  first  attempt  is  made  to  use  seismic  data  in  gravity  prediction.  The  out¬ 
lined  approaches  are  rather  possible  proposals  than  fully  developed  algorithms.  Still  more 
work  has  to  be  done  in  that  field  since  geophysicists  started  to  develop  three-dimensional 
digital  seismic  velocity  models  of  (parts  of)  the  earth  which  will  be  available  in  near  future. 
In  our  study  a  corresponding  model  was  developed  for  the  European  Alps  (see  chap.  8). 


7.1  Simple  seismic  velocity-density  relationship  in  a  homogeneous  medium 

Through  seismic  experiments  artificial  shocks  (or  explosions)  are  generated  to  which  the 
earth  responds  by  releasing  part  of  the  energy  in  form  of  elastic  waves  which  travel  through 
rocks  with  a  certain  velocity  depending  on  density  and  elastic  moduli.  The  same  happens 
when  earthquakes  occur.  The  types  of  body  waves  (there  are  also  surface  waves)  follow  the 
laws  of  geometrical  optics  being  reflected  and  refracted  at  layer  boundaries  where  the 
velocity  (and  consequently,  the  material  constants)  changes. 

Assuming  a  homogeneous  medium  we  observe  the  so-called  P  waves,  longitudinal  waves,  due 
to  transmission  of  compressions  and  rarefactions,  whose  velocities  vp  are  given  by 


where  A,  p  are  the  Lame  elastic  parameters  characterizing  the  material,  and  p  is  the  density. 
In  seismology,  it  is  convenient  to  replace  the  parameter  A  by  the  so-called  bulk  modulus  k, 

k  =  A  +  |m  (7-2) 


so  that  we  get  for  P  waves 

k  +  ^l0-5 

~p~ 

li  is  also  called  the  rigidity  modulus. 


(7-3) 
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The  S  (shear)  wave  velocities  (transverse  wave  motion)  are  given  by 

vs  =  (p/p)0-5-  (7_4) 

Thus,  (7-3)  and  (7-4)  represent  the  relationship  between  seismic  velocities  and  density,  so 
far  the  bulk  and  rigidity  moduli  are  known  and  a  homogeneous  medium  is  assumed.  The 
non-linear  relations  can  be  linearized  and  -  at  least  theoretically  -  observation  equations 
can  be  built  up  as  a  function  of  the  disturbing  gravity  potential  using  Poisson’s  equation 
(6-1).  Remember,  however,  the  implications  with  it  discussed  in  chap.  6.  By  doing  so,  we  get 
(A  is  here  the  Laplace  operator) 

<5vp  =  (8 *G)-‘  (k0  +  Imo)0'5  Pa''  “  AT  (7-5) 

and 

<5vs=  (8irG)-' p0p~oUS  AT  (7-6) 

where  AT  indicates  that  the  relations  to  second  order  derivatives  AT  =  £  dT/dxi  are 

1=1 

established.  tfv„  =  vp  -  vp  ,  6vs  =  vs  -  vs  are  the  residual  observations,  and  k0,  po,  po  are 

o  o 

known  (approximate)  values. 

As  extensively  discussed  in  chap.  3  and  verified  by  own  calculations,  we  can  find  a  linear 
empirical  relationship  between  densities  and  P  wave  velocities  of  the  form  (BIRCH  1961) 

p  =  a  +  S  Vp  (7-7) 

or,  expressed  as  vp  =  f(p)  ,  we  get 

v,,  =  -a  +  bp  (7-8) 

where 

b  =  b'1  =  dv/dp,  a  =  ba  (7-9) 

BIRCH  (1961)  has  published  for  a  =  0.41  [g  cm'3}  and  for  6  =  0.3597  [(g  cnr3)/(km  sec'1)] 
corresponding  to  approximately  b  =  2.8  [(km  sec'>)/(g  cm"3)]  which  he  found  from  laborato¬ 
ry  experiments.  It  is  well  known  from  studies  with  actual  data,  that  those  coefficients  b 
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found  are  somewhat  higher  in  the  range  4  ...  5.5  [(km  sec1)/(g  cm'3)].  Our  own  experiment 
showed  a  mean  value  of  5.4  ±  0.4  for  the  European  Alps.  Even  more  complex  relationships 
can  be  established;  see,  e.g.,  ANDERSON  (1967,  1970),  BUNTEBARTH  (1982),  RY- 
BACH,  BUNTEBARTH  (1982). 

Thus,  assuming  that  a  sufficient  trend  (or  constant)  with  respect  to  the  vp  velocities  is 
subtracted  from  the  data,  so  that  we  can  work  with  residual  or  anomalous  velocities  ftv 
forming  a  scalar  stationary  random  process,  the  coefficient  b  (7-9)  can  be  used  for 
constructing  the  relevant  covariance  propagation  which  is  needed  to  consider  it  in  the 
integrated  model  for  gravity  prediction. 

We  further  think  that  digital  three-dimensional  seismic  velocity  models  will  be  available  in 
near  future  -  we  have  one  constructed  over  the  European  Alps,  see  chap.  8  -  so  that  those 
informations  can  be  used  to  determine  and/or  to  improve  the  gravity  field.  In  this  context, 
also  proposals  and  developments  for  global  P  velocity  models  in  the  form  of  spherical  har¬ 
monic  models  are  already  made  (DZIEWONSKI  1984;  PEC,  MARTINEC  1985;  WOOD- 
HOUSE,  DZIEWONSKI  1984). 

7.2  Seismic  wave  motion  in  an  inhomogeneous  medium 

Since  density  p  and  elastic  Lame  parameters  A,  p  (or  bulk  and  rigidity  moduli,  respectively) 
are  very  irregular  junctions  of  position,  they  may  be  treated  as  randomly  distributed  in 
space  forming  a  stochastic  process.  These  statements  refer  to  the  famous  work  of  CHER¬ 
NOV  (1960).  Meanwhile  among  geophysicists  the  theory  was  successfully  applied,  see,  e  g. 
AKl  (197S),  CAPON  (1974),  KNOPOFF,  HUDSON  (1964),  KORN  (1987)  -  to  mention 
also  some. 

Thus,  KNOPOFF,  HUDSON  (1964)  introduced  six  correlation  functions  of  the  form 


Cxx 

=  cov(«5A(P),  <5A(Q)) 

(7-10) 

Cp, 

=  cov(<%P),  6p(Q)) 

(7-11) 

Cpp 

=  cov(6p(P ),  6p( Q)) 

(7-12) 

CXp 

-  cov(<iA(P),  <V(Q)) 

(7-13) 
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'Xp 


»4> 


=  cov(£A(P),  £p(Q)) 
=  cov(^P),  6p(Q)) 


(7-14) 

(7-15) 


and  treated  the  quantities  6p,  6X,  6p  as  signals  which  perfectly  fit  to  our  proposed  integra¬ 
ted  approach  in  conformity  to  least-squares  collocation. 

Following  KNOPOFF,  HUDSON  (1964)  the  basic  equation  of  wave  motion  in  an  inhomoge¬ 
neous  elastic  medium  is  given  by  the  following  partial  differential  equation 


d2u 


+  2  (V  p  •  V)  h  -  2  (V  p)  V  •  u  +  2  (V  p)  *  V  *  u 


(7-16) 


where 

a 

P 

P,  X 

t 

V 

V- 

V* 


seismic  displacement  vector  (observed  by  seismographs), 
density, 

elastic  parameters  (Lame  constants), 
time, 

gradient  operator, 
divergence  operator, 
rotation  operator. 


For  the  sake  of  simplicity  we  want  to  discuss  only  the  effect  of  density  anomalies  on  wave 
motion.  Thus,  we  assume,  that  X  and  p  are  (known)  constants.  This  assumption  does  not  al¬ 
ways  conform  with  reality,  but  simplifies  the  discussion.  We  are  mainly  interested  in  densi¬ 
ty  anomalies,  since  their  relationship  to  gravity  was  discussed  in  the  chapters  before. 

The  generalization  with  respect  to  all  three  parameters  p, X,p  can  be  done  using  the  work  of 
KNOPOFF,  HUDSON  (1964)  and  CHERNOV  (1960).  The  formal  solution  already  exists. 

Starting  point  is  again  the  decomposition 


p  =  p0+  6p 

X  =  Xq 
P  =  Po 


(7-17) 

(7-18) 

(7-19) 
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where  p0,  X0,  Po  are  approximate  values  and  constants,  and  6p  are  density  anomalies.  In  the 
same  way  we  decompose  the  displacement  vector  u  in 

u  =  u<)  +  £u  (7-20) 


where 

uq  is  the  approximate  value,  and 

<5u  is  the  variation  of  the  displacement  vector  due  to  density  anomalies. 

As  a  consequence  of  (7-17)  to  (7-19)  the  gradients  are  VAo  =  Q,  Vpo  =  Q,  and  =  Q,  which 
enter  the  basic  equation  of  motion  (7-16).  It  can  be  linearized  resulting  in 

-I.)-® *■(»■«.) 

mi  = 

where 

f  =~^6P  (7‘23) 

(7-21)  defines  the  approximate  wave  motion,  if  the  medium  would  be  homogeneous. (7-22) 
is  the  inhomogeneous  elastic  wave  equation  which  relates  the  variation  tfu  of  the  displace¬ 
ment  vector  to  the  density  anomalies  6p,  where  f  is  the  disturbing  force. 

It  is  further  common  to  introduce  the  wave  velocities  a  and  0  defined  by 

a2  =  (7-24) 

P  o  v  ' 

0i  =  U*  (7-25) 

P  o 

For  the  solution  of  wave  equations  like  (7-21),  (7-22)  we  decompose  the  displacement  vec¬ 
tor  in  one  part,  which  is  due  to  a  compression  (P)  wave,  and  in  another  part,  which  is  due 
to  a  shear  (S)  wave.  Thus,  introducing 

6<(>  =  V  •  <5u  (7-26) 

Sw  =  V  »  <5u  (7-27) 


(7-21) 

(7-22) 
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we  find  from  (7-22) 


d36v 


=  a3  A  6<f>  +  V  •  i 

P  wave 

(7-28) 

=  (31  A  6w  +  7  *  f 

S  wave 

(7-29) 

where  A  is  the  Laplace  operator. 

Assuming  we  have  solved  (7-28),  (7-29)  the  diplacement  vector  variation  6u  is  given  by 

5u  =  V  6<f>  +  V  *  £w  (7-30) 

Under  the  assumption  that  the  incident  WLve  uo  is  harmonic,  that  means  time-dependent 
with  e‘iu,t,  the  solutions  for  6<p  and  Sw  can  be  found  in  KNOPOFF,  HUDSON  (1964). 


{*=-nrb?/(,'-£)Lr-dv 

V 

(7-31) 

(7-32) 

V 


Where  (r,  0,  A)  are  spherical  coordinates  (radius,  co-latitude,  longitude)  and  ^  is  the 
spherical  distance  (P,Q) 


1  = 

(r*  +  i' 3 

-  2  r  r'  cos^)]°- 5 ,  r  =  r(P),  r'  =  r(Q) 

(7-33) 

dv  = 

t' J  sin  O' 

dr' d  0’  dA 

(7-34) 

II 

e 

w/ a 

(7-35) 

k0  - 

“IP 

(7-36) 

V'  = 

3Fer  + 

1  d  _  ,  d  _ 

r'  sin 0'  di'  *  +  r'  ffl 

(7-37) 

(gradient  operator  in  spherical  coordinates) 

w  is  the  frequency  of  the  incident  wave,  and  using  (7-31),  (7-32)  we  find  for  the 
displacement  vector  variation  fiu  (7-30) 


+  nrwv’f  (fS-r-)  ■(''  ■0dv  V-W 
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We  notice  the  important  fact,  that  the  integral  solution  for  (7-38)  is  in  principle  of  the 
same  form  as  the  expressions  for  the  gravity  field  functionals  using  Newton’s  attraction  in¬ 
tegral.  The  variation  <$u  of  the  displacement  vector  is  also  given  as  an  integral  over  the' 
volume  v,  where  v  should  coincide  with  the  volume  in  which  we  have  density  anomaly  data. 
The  density  anomalies  enter  in  (7-38)  via  the  disturbing  vector  f. 

The  incident  wave  uo  is  assumed  to  be  harmonic,  e  g.,  Uo  is  a  function  of  the  form 


u0(P,t)  =  e‘iu,t  u0(P) 


(7-39) 


with  time-dependence  e_ililt  and  P(r,  0,  A).  Considering  (7-39)  we  find  for  f  (7-23) 

f(P,t)  =  w2  e*iu>t  5p(P)  Uo(P)  .  (7-10) 

In  order  to  evaluate  the  integrals  for  the  displacement  vector  bu  we  have  to  insert  f  at  the 
integration  point  Q(r',  O' ,  A')  into  the  integrals  and  to  compute  the  divergence  and  rota¬ 
tion.  Thus, 

(V'  •  f)  =  [V'  Sp( Q)  u0(Q)  +  6p( Q)  V'-  u„(Q)]  (7^1) 

(V'  *  f)  =  v2e- i*t  (V'  6p( Q)  »  uo(Q)  +  Sp{ Q)  V'*  u0(Q)]  (7-42) 

We  can  see  from  (7-41),  (7-42)  that  the  displacement  vector  variation  bu  is  not  only  a  Junc¬ 
tion  of  bp  itself,  but  also  a  function  of  the  density  gradient  V'  bp  .  This  means  that  a  change 
of  density  in  the  volume  v  will  cause  an  effect  on  the  displacement  vector,  a  physically  real¬ 
istic  model. 


For  the  components  of  V 
spherical  coordinates: 
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,  we  get  the  following  expressions  using 

(7-43) 

(7-44) 

(7-45) 
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The  derivatives  of  (1/1)  are 


(7-46) 

(7-47) 

(7-48) 


7.3  A  possibility  of  evaluating  the  autocovariance  function  of  density  anomalies  from 
seismic  observations 

The  concept  of  solution  assumes  that  we  know  the  autocovariance  function  of  density 
anomalies  C6p6p.  The  computation  can  be  done  using  available  density  data  as  outlined  in 
chap.  6  and  done  numerically  in  this  study  (see  chap.  8,  9). 

Another  method  of  determination  of  Cgp6p  =  cov(6p,  6p)  is  outlined  in  (KNOPOFF,  HUD 
SON  1964)  which  is  similar  to  variance-component  estimation.  The  method  is  based  on  the 
analysis  of  seismic  observations  1  =  (£u(P,t))  providing  an  a  priori  base  function  repre¬ 
sentation  of  C6p6p  with  unknown  com,'  mts.  KNOPOFF,  HUDSON  (1964)  have  used  a 
Gauss’  function  in  their  approach. 

C6p6p(p>Q)  =  Cgpe~aP  IP-QI  (7-57) 

But  we  may  also  think  of  a  spherical  harmonic  representation  of  C6p6p  : 

N 

CWP>Q)  =  X  Cn(p>Q)  ‘  Pn(cos^  (7“58) 

n  =  0 

In  (7-58)  the  factors  c„(P,Q)  may  be  treated  as  constants  (in  very  local  investigations); 
they  may  be  some  functions  of  radial  distances,  considering  the  variations  of  density  with 
depth. 

Let  us  assume  now,  that  we  have  a  network  of  seismic  stations  Pi  at  the  earth’s  surface.  At 
each  seismic  station  Pi,  it  is  possible  to  observe  the  displacement  vector  1  =  (£u(Pi,t))  as  a 
function  of  time  for  a  seismic  event. 

With  1  we  start  the  covariance  analysis  by  computing  for  each  station  Pi  the  covariance 
matrix  Cn  =  (Pi,Pi) 


Cn  =  E  {  11*}  (7-59) 

Note,  that  *  denotes  here  the  conjugate-complex  transposition,  because  we  deal  with 
complex  operators  and  observations. 
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Using  (7-63)  in  (7-59)  we  get 


fiu(Pi,Pi)  =  AG*p8pA*  +  A  1*  +  B  A*  +  A  B*  +  Cn„  (7-60) 

Having  assumed  that  the  incident  wave  Uo  is  harmonic  with  time-dependence  e*iu)t  we 

notice,  that  for  a  station  Pj  1  1*  is  now  no  more  any  longer  time-dependent  because  of 
£*ilDt  .  gitfit  —  X. 

The  computation  of  £n  is  obtained  by  evaluation  of  the  seismograms  in  vertical  and 
horizontal  direction. 

The  matrix  £jp8p  consists  of  all  covariances  between  the  discrete  points  in  the  volume  v. 
The  further  way  is  to  substitute  a  model  for  a  covariance  function  like  (7-57)  or  (7-58) 
into  the  covariance  matrix  C^p.  For  example,  the  covariance  function  (7-57)  is  a  function 
of  the  distance  |  P— Q  |  with  two  free  parameters,  the  variance  Cgp  and  the  distance-scale 
parameters  ap  to  be  determined. 

Under  the  assumption,  that  we  know  some  approximate  values  for  C8p,  ap,  we  may  expand 
Cn  ,  see  (7-60),  into  a  Taylor  series,  neglecting  all  higher-order  terms, 

£n(P|.  Cf,,  op)  =  fil, (Pi,  C(»p,  <$  +  ^  K  (7-61) 

Using  the  observation  vector  1,  we  may  compute  the  so-called  mean  square  fluctuation 
(KNOPOFF,  HUDSON  1964),  that  is  the  variance  Co(l). 

C0(l)  =  f  £n  1  (7-62) 

Considering  (7-61)  together  with  (7-62)  it  should  be  possible  to  minimize  the  difference 
SCt(Pi) 

6Ci(Pi)  =  C0(l)  -  Cj(l)  =  1*  Cn  1  -  f  Cu  l  +  f  ^  l  *Cgp  +  1*  ^  1  6ap  (7-63) 

with  respect  to  0C£p  and  6ap  (the  unknown  parameters  of  the  covariance  function  (7-57) 
for  all  seismic  stations  Pi,  i  =  {1,  ...  ,  N}. 
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In  this  way  we  may  obtain  improved  parameters  for  the  autocovariance  function  Cbpi)p 
Eqs.  (7-59)  to  (7-63)  are  in  principle  a  discrete  approach  to  the  discussions  of  KNOPOFF, 
HUDSON  (1964).  Here  an  integral  expression  for  the  mean  square  fluctuation  (7-62)  is  de¬ 
rived.  Subsequently,  the  covariance  function  (7-57)  is  inserted  into  the  integral  and  the 
integration  is  performed.  The  result  of  integration  yields  also  the  mean-square  fluctuation 
as  a  function  of  Cgp  and  ap  . 
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8.  TEST  DATA 


8.1  Test  network  "Rossdorf 

The  test  area  Rossdorf  is  located  around  50  km  south  of  Frankfurt /F.R.  Germany.  It  has  an 
extension  of  26  *  23  km2  and  is  especially  well-suited  because  of  large  variations  in  (surface) 
density.  This  is  due  to  the  fact  that  part  of  the  area  belongs  to  the  Rhinegraben  whose 
upper  layer  consists  of  sedimentary  fillings.  The  eastern  part  of  the  test  area  belongs  to  the 
so-called  "Odenwald",  a  low  mountain  range.  For  the  investigations  78  (point)  gravity  ob¬ 
servations  were  available.  At  these  stations,  at  least  in  the  near  neighbourhood,  also  surface 
densities  were  determined  by  a  geologist  (FAHLBUSCH  1987,  pers.  comm  ).  The  density 
data  are  estimates  coming  from  laboratory  investigations,  boreholes,  and  descriptions  found 
in  corresponding  archives.  They  should  be  representative  for  a  cylinder  of  about  150  m  in 
diameter  around  the  gravity  station  and  a  depth  of  approximately  1  km. 

The  coverage  of  the  area  with  observational  points  cannot  be  described  to  be  homogeneous, 
because  they  follow  more  or  less  levelling  lines,  see  Fig.  8.1.  The  density  values  vary  largely 
from  2,48  g/cm3  to  3.02  g/cm3.  The  area  has  only  height  variations  of  a  few  hundred 
meters.  The  topography  is  plotted  in  Fig.  8.2  and  a  contour-line  plot  of  surface  densities  is 
given  in  Fig.  8.3. 
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8,2  Test  area  "European  Alps" 


The  test  area  European  Alps  extends  from  45.0°  to  48.5°  North  and  from  6.0°  to  16.5°  East. 
It  comprises  Switzerland,  most  of  Austria,  southern  Germany  to  the  Danube  River, 
northern  Italy  including  the  Po  plain,  and  easternmost  France  excluding  the  western  Alps 
and  the  Lake  Alps.  The  following  section  describes  the  data  that  are  used  in  the  geophysical 
modelling  (see  chap.  3)  as  well  as  in  the  numerical  computations  with  respect  to  gravity 
prediction.  Part  of  the  data  are  digitally  represented  on  an  equidistant  grid  of  6'  *  10' 
(Ay?  *  AA),  corresponding  to  an  area  of  11  *  13  km2.  The  grid  was  chosen  because  other 
data  sets  were  already  available  on  it  (area  means  of  topography  and  free-air  anomalies). 
The  data  values  are  considered  as  mean  values  for  the  grid  compartments.  Gridded  data 
used  in  our  investigations  are  digital  models  of  topography  (digital  terrain  models), 
equiangularly-spaced  mean  block  values  of  gravity  anomalies,  and  a  model  of  surface  rock 
densities. 

Block  mean  values  of  gravity  and  point  values  have  considerably  different  spectral  proper¬ 
ties.  There  was  also  a  need  to  use  point  values  of  gravity  field  functionals,  especially  in  the 
context  of  least-squares  collocation  for  the  determination  of  a  corresponding  reasonable 
covariance  function  of  the  gravity  anomalies. 


6.2.1  Topography 

The  anomalous  gravity  field  is  primarily  due  to  two  facts  that  seem  at  the  first  glance  quite 
different:  variations  in  the  visible  topography  and  radial  as  well  as  lateral  inhomogeneities 
of  the  sub-surface  masses.  A  major  part  in  local  gravity  field  variations  (medium  and  high 
frequencies  of  the  spectrum  of  the  gravity  field)  is  a  direct  consequence  of  the  topography. 

So  one  of  the  fundamentals  of  a  proper  modelling  of  the  gravity  field  is  a  model  of  the 
heights  in  the  area  of  interest.  Heights  stored  in  gridded  form  together  with  the  assumption 
of  a  constant  density  form  a  zero-order  model  of  the  distribution  of  masses  whose  successful 
use  has  been  demonstrated  in  various  applications  (e.g.,  KEARSLEY  et  al  1985;  FORS- 
BSRG,  TSCHERNING  1981). 

For  our  computations  several  sets  of  height  data  were  collected  and  used: 

-  The  global  set  of  block  mean  values  of  heights/depths  released  by  the  working  group  of 
Prof.  Siinkel,  Technical  University  of  Graz,  described  in  (WIESER  1988).  This  set, 
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referred  to  as  TUG87m5,  has  a  resolution  of  5'  by  5'  and  is  referring  to  geographical  co¬ 
ordinates  (<£>,  A).  It  is  compiled  from  two  sets  of  5'  by  5'  resolution:  ET0P05  and 
DB5B5  (previously  SYNBAPS).  Processing  of  both  sets  consisted  not  only  of  merging 
the  data  but  primarily  included  also  detection  and  removal  of  gross  errors. 

The  Austrian  digital  terrain  model  with  a  resolution  of  20"  by  20"  is  based  on  geographi¬ 
cal  coordinates  and  was  developed  by  the  Institute  of  Photogrammetry  of  the  Technical 
University  of  Vienna  (HAITZMANN  1983).  This  model  is  referred  to  as  POHEKR,  see 
Fig.  8.4  for  a  plot  of  a  1°  *  1°  block. 

For  Switzerland  a  set  of  approximately  100  m  by  100  m  point  values  was  available  called 
RIMINI  ( De  MARCHI  1983).  The  model  is  based  on  the  map  1  :  25  000  of  the  Federal 
Office  of  Topography.  It  comprises  whole  Switzerland  (a  total  of  330  maps  1  :  25  000) 
resulting  in  1.15  million  height  values. 
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Fig  8.4.  Heights  in  a  1°  «  1°  block  of  the  European  Alps. 
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->  For  the  southern  part  of  Germany  a  set  of  point  values  based  on  UTM-coordinates  is 
available.  Its  resolution  is  approximately  100  *  100  m2.  For  a  limited  area  we  had  also 
access  to  a  model  of  resolution  -  50  *  50  m2. 

An  overview  over  all  available  sets  is  given  in  Fig.  8.5. 


8.2.2  Gravimetric  observations 

The  largest  data  set  of  observations  consists  of  an  amount  of  41  681  gravity  values  within 
our  Alpine  area.  For  an  overview  the  data  have  been  thinned  to  a  minimum  point  separa¬ 
tion  of  10  km  and  plotted,  see  Fig.  8.6. 

The  gravity  values  for  the  eastern  part  of  France  as  well  as  for  whole  Switzerland  have  been 
made  available  by  the  Bureau  Gravimetrique  International,  Toulouse.  The  Italian  petrol 
company  AG  IP  provided  a  data  set  covering  that  part  of  northern  Italy  which  was  of 
interest  for  our  investigations.  The  Technical  University  of  Graz,  Austria,  provided  the 
data  for  the  Austrian  part  of  our  area. 

The  computation  of  the  mean  Bouguer  anomaly  values  for  the  6'  *  10'  areas  is  based  on 
the  mean  free-air  anomaly  and  elevation  values  for  Europe,  as  supplied  by  Institut  fur  An- 
gewandte  Geodasie  (IFAG).  Application  of  the  mean  Bouguer  reduction  to  the  mean  free- 
alr  anomalies  results  in  an  anomaly  which  approaches  the  mean  Bouguer  anomaly.  In  the 
Alps  this  anomaly  is  strongly  correlated  (with  negative  sign)  with  the  mean  elevation.  It  re¬ 
flects  the  fact  that  the  mean  elevation  is  closely  correlated  with  Moho  depth.  The  simple 
Bouguer  reduction  applied  is  defined  by: 

6gn  =  2ir  G  p  H  (8-1) 

with  H  being  the  mean  elevation  value  of  a  compartment,  G  =  6.673  •  10‘8  [g‘>  cm3  s'2], 
5gu  in  mGal.  The  density  p  was  everywhere  assumed  to  be  2.67  g/cm3,  the  value  usually 
taken  for  the  average  crust. 

For  relatively  flat  areas  the  approximation  of  an  infinitely  extended  Bouguer-plate  is 
sufficient  but  for  most  of  our  area  the  full  terrain  reduction  had  to  be  computed.  This  was 
done  using  the  prism  integration  method  which  allowed  also  to  simultaneously  determine 
the  Isostatic  effect  due  to  the  topographic  load.  For  the  isostatic  compensation  the  Airy- 


78 


MUSE 


T 


Fig.  8.6.  Available  point  gravity  anomalies  over  the  European  Alps. 
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Heiskanen  compensation  mechanism  was  used  assuming  the  following  parameters:  p  =  2.67 
g/cm3,  T  =  32  km,  Ap  =  0.4  g/cm3. 

Unfortunate!/  there  are  gaps  in  the  given  data  sets  resulting  in  significant  gaps  of  the 
Bouguer  anomaly  particularly  in  the  Central  Alps.  Fig.  8.7  is  a  crude  representation  of  the 
distribution  of  Bouguer  anomaly  values. 


8  2.3  Astronomical  observations  (deflections  of  the  vertical) 

For  Austria  a  set  of  astronomical  latitudes  and  logitudes  (deflections  of  the  vertical)  with  a 
total  of  521  values  was  available.  Field  work  and  accuracies  are  documented  in  ( BRAND - 
ST  A  TTER  1987).  Parts  of  the  observations  were  carried  out  by  the  Institute  of  Theoretical 
Geodesy  and  Geophysics,  Vienna,  Austria  (Prof.  Bretterbauer). 

Switzerland  (Prof.  Kahle,  ETH  Zurich)  provided  a  set  of  240  deflections  of  the  vertical.  A 
good  number  of  them  are  published  in  (GURTNER  1978). 

For  the  southern  part  of  Germany  the  deflections  of  the  vertical  are  taken  from 
( SCHMIDT ,  EHLERT  1982). 

For  a  plot  of  all  available  deflections  of  the  vertical  see  Fig.  8.8. 


8.2.4  Seismic  data 

The  crustal  structure  of  the  Alpine  region  has  been  thoroughly  investigated  mainly  by 
seismic  refraction  experiments  during  the  last  30  years.  The  results  have  been  published  in 
a  large  number  of  papers  and  they  do  not  cover  the  body  of  the  Alps  completely.  In  recent 
years  a  number  of  attempts  have  been  made  to  collect  all  the  data  and  to  arrive  at  a  unified 
areal  interpretation  (e  g.,  MOSTAANPOUR  1984;  GEISS  1987).  In  spite  of  the  doubtless 
merit  of  these  attempts,  we  have  not  simply  incorporated  the  results,  partly  because  we  are 
aiming  at  greater  detail,  partly  because  the  uncritical  use  of  other  workers’  results  is  easily 
misleading  For  our  computations  of  the  gravity  and  geoidal  effects  of  the  seismically 
inferred  crustal  structures  we  have  again  reviewed  and  analysed  all  the  original  publications 
available  to  us. 
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The  crust  and  upper  mantle  have  a  complicated  three-dimensional  structure  which  cannot 
possibly  be  represented  adequately  except  on  a  very  fine  three-dimensional  grid.  However, 
there  are  elements  of  order  to  which  the  various  seismic  techniques  are  sensitive.  Refraction 
experiments  (including  the  methods  of  data  analysis  and  interpretation)  will  tend  to  "see" 
refractors  or  layers,  while  reflection  experiments  mainly  "see"  reflectors  or  surfaces  of 
seismic  impedance  contrast  mostly  on  a  much  smaller  scale  (vertically  and  horizontally); 
some  may  coincide  with  refractors,  but  most  of  them  reflect  internal  structures  of  the  gross 
layering  derived  refraction  seismology. 

For  convenience  and  as  a  means  of  data  condensation  we  divided  the  crust  into  two  layers 
and  the  upper  mantle  similarly  into  the  lithospheric  part  and  the  underlying  asthenosphere 
(with  no  bottom  assumed).  The  most  important  refractor  or  large-scale  boundary  is  the 
Mohorovicic  discontinuity  dividing  crust  and  upper  mantle;  its  depth,  or  equivalently  the 
crustal  thickness,  Zm,  is  reasonably  well  established  in  the  region  of  study.  Another 
boundary  on  which  we  have  significant  information  is  the  top  of  the  crystalline  basement  or 
its  depth,  Zb-  The  depth  of  the  lithosphere-asthenosphere  transition  (or  "boundary"),  Zj,  is 
known  with  much  less  certainty;  it  is  based  mainly  on  the  dispersion  of  (earthquake)  surface 
waves  with  poor  spatial  resolution. 

To  start  with,  our  seismic  model  can,  thus,  be  represented  by  three  layer  boundaries  at 
depths  Zb,  Zm,  and  Zi.  However,  as  argued  above,  the  layers  are  not  homogeneous;  in  the 
individual  refraction  results  there  is  ample  evidence  for  seismic  velocity  variations  within 
the  layers.  We  attempted  to  treat  these  (as  far  as  the  data  permitted  it)  as  a  perturbation 
of  the  layered  model. 

First  we  constructed  the  layer  boundaries  on  the  basis  of  published  contour  maps  and, 
particularly,  of  seismic  sections  along  observed  profiles  and  a  few  fans.  The  sections  usually 
show  boundaries  and  seismic  P  wave  velocities,  vp,  mainly  at  the  top  of  the  layers.  The 
refraction  results  are  also  often  published  as  velocity-depth  functions  for  refraction  profiles 
or  parts  of  them,  pertaining  to  certain  regions.  In  addition  reflection  seismic  results  have 
been  used.  As  mentioned,  the  lithosphere  bottom  was  taken  from  interpretations  of  surface 
wave  dispersion  and  partly  of  teleseismic  traveltime  residuals 

All  these  data  were  digitized.  The  irregular  points  were  then  interpolated  onto  the  grid 
described  above.  Digitizing  the  input  data  had  to  take  into  account  the  different 
cartographic  projections,  so  that  it  was  done  mostly  manually  with  the  aid  of  strongly 
magnified  maps.  To  reduce  the  errors  of  locations  (of  the  data  points),  mostly  points  of 
intersection  of  contours  with  lines  of  latitude  or  longitude  were  taken.  At  the  beginning, 
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only  those  results  were  accepted  that  had  been  considered  reliably  by  the  original  authors, 
and  extrapolations  to  the  marginal  regions  (mostly  shown  by  dashed  lines)  were  omitted. 
Only  in  larger  regions  with  no  data  we  included  such  information  of  less  reliability.  In 
regions  of  data  overlap  from  two  or  more  sources  we  tried  to  avoid  discrepancies  in  the  first 
step.  In  the  second  step  the  data  set  was  analysed  with  respect  to  significant  discrepancies. 
These  were  either  eliminated  altogether  or,  in  regions  of  broad  scatter,  the  arithmetic  mean 
was  taken. 

This  procedure  resulted  in  a  set  of  irregularly  distributed  points  values.  In  the  third  step 
this  data  set  was  interpolated  onto  a  regular  grid.  Two  methodes  were  used:  the  method 
due  to  Hardy  (GOPFERT  1977),  called  "multiquadratic",  has  been  applied  if  not  otherwise 
stated;  alternatively  interpolation  by  triangulation  was  used  (AKIMA  1978),  were  a  5th 
order  polynomial  is  expanded  within  triangles.  The  multiquadratic  interpolation  renders  a 
smoother  set  of  point  values  and  it  is  less  sensitive  to  greater  data  gaps  and  strong 
gradients  resulting  from  discrepancies  between  close-by  points;  however,  computation  takes 
much  longer  so  that  we  were  limited  in  the  numbers  of  input  and  output  points  because  of 
limitations  enforced  by  the  reasonable  computing  time.  The  resulting  grid  (75-120  %  of  the 
compartment  size,  depending  on  the  original  data  density)  was  further  refined  to  20  %  of 
the  compartment  size  with  the  aid  of  bicubic  splines  (PRESS  et  al  1986).  Subsequently  all 
grid  points  thus  obtained  and  lying  within  a  6'  *  10'  quadrangle  were  arithmetically 
averaged  and  the  mean  value  is  then  taken  for  this  quadrangle.  Averaging  naturally  results 
in  some  smoothing  of  the  data,  but  methodical  weakness  of  the  different  kinds  of  interpola¬ 
tion  is  partly  cancelled.  Tests  in  which  both  interpolation  methods  were  compared  for  some 
regions  have  demonstrated,  that  the  results  after  the  final  averaging  are  nearly  identical. 
Thus,  we  believe  that  the  requirement  of  a  uniform  data  generation  has  been  satisfied. 

Lateral  variations  or  perturbations  of  the  seismic  velocities  within  the  layers  were  estimated 
for  parts  of  the  crust  in  the  following  manner.  We  started  from  published  velocity-depth 
functions  v(z)  given  for  certain  "points".  Assuming  a  piecewise-linear  velocity  distribution, 
we  computed  mean  velocity  v,  for  n  layers  i  of  thickness  The  average  velocity  v*  of  a 
crustal  section  Z  =  E  Zj  is  given  by  (MOSTAANPOUR  1984): 

n 


S  Zj 


Fig  8.9  shows  the  location  of  crustal  sections  discussed  beside  the  contour  maps.  In  the 
following  the  seismic  boundaries  are  presented  from  top  to  bottom,  followed  by  the  lateral 
velocity  perturbations. 
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8.2.5  Depth  of  the  seismic  basement 


The  seismic  basement  is  defined  by  the  depth  Zi,  where  the  P  wave  velocity  vp  first  reaches 
the  value  6.0  km/s,  not  by  reflecting  horizons  as  usually  done  in  exploration  geophysics.  If 
low-velocity  layers  exist  at  greater  depth  there  may  be  additional  6.0  km/s  surfaces  which 
are  not  included  in  the  "seismic  basement".  What  is  the  significance  of  the  basement  as  de 
fined  here7  In  a  general  way,  it  will  correspond  to  the  top  of  the  crystalline  rocks;  thus  the 
depth  Z|,  is  equivalent  to  the  thickness  of  the  sedimentary  "layers".  This,  however  is  not 
everywhere  to  be  taken  literally,  particularly  not  in  the  Alps.  Shallow  crystalline  rocks,  par- 
ticularly  if  tectonically  deformed  usually  have  vp-values  below  6.0  km/s.  This  has  been  de¬ 
monstrated  by  BIRCH  (1958)  for  normal  granite  assuming  a  temperature  gradient  of  18  °C 
per  km  (Fig.  8.10).  The  relatively  low  velocities  are  caused  by  porosity  and  particularly  by 
small  open  joints  related  to  weathering  and  tectonics;  at  greater  depth  pressure  generally 
closes  the  joints  and  vp  exceeds  6.0  km/s,  usually  at  depth  of  about  2  km.  There  are, 
however,  also  sound-hard  sediments,  e.g.,  dolomite,  with  vp-values  greater  than  6.0  km/s. 

The  seismic  basement  Z|,  as  assumed  here  was  essentially  based  on  the  interpretation  and 
the  contour  maps  published  by  MOSTAANPOUR  (1984)-  In  addition,  we  took  a  number  of 
velocity-depth  functions  at  a  few  points  along  the  profiles  shown  in  Fig.  8.9;  these  points 
were  incorporated  into  the  interpolation  or  were  used  for  checking.  The  result  of  our 
interpolation  is  presented  in  Fig.  8.11. 

The  contours  reflect  the  tectonic  gross  structures  of  the  study  area.  The  evolution  of  the 
Alps  that  resulted  in  these  structures  may  be  briefly  characterized  as  follows.  A  number  of 
compressive  pulses  lead  to  folding  and  nappe  movement  from  the  central  region  towards  the 
northern  foreland.  Today  the  Alps  have  an  asymmetric  structure  of  the  deformed  crust. 
While  in  the  north  one  crustal  nappe  after  the  other  separated  from  its  basement  and 
moved  north,  in  the  south  a  suture  formed,  the  Insubric  line.  It  separates  the  East  Alpine 
Block  from  the  South  Alpine  block.  To  the  north  and  south  there  are  two  large  sedimentary 
basins,  the  Molasse  and  the  Po  basin,  which  received  the  erosional  products  of  the  Alpine 
orogene  during  repeated  uplift  phases.  In  east-west  direction  the  Helveticum  and  the 
Penninicum  are  to  be  distinguished.  Fig.  8.12  shows  the  tectonic  gross  structures  of  the 
Swiss  Alps. 

The  greatest  values  of  Zb  are  formed  in  the  sedimentary  basins  of  the  Molasse,  the  southern 
Rhine  Graben,  and  the  Po  region.  In  the  Molasse,  the  maximum  values  of  7  to  10  km  are 
reached  at  the  edge  of  the  Alps.  In  the  southern  Rhine  Graben,  values  up  to  11  km  are 
reached  between  the  much  shallower  values  in  Black  Forest  and  Vosges  massifs.  The  Po 
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basin  shows  up  with  values  up  to  8  km.  The  shallowest  depths  of  the  basement  (less  than 
2  km)  are  found  in  the  central  East  Alpine  regions  and  along  a  narrow  arch  of  the  western 
Alps  (the  zigzag  shape  of  the  contours  results  from  insufficient  grid  resolution).  Somewhat 
greater  Zb  values  (3  ...  5  km)  correspond  to  the  Helveticum  and  Penninicum,  rising  towards 
the  Molasse  basin.  In  the  South  Alpine  region  the  Zb  values  are  variable  between  2  and 
4  km  before  they  drop  towards  the  sedimentary  Po  basin.  The  transition  is  here  more 
gentle. 


5,2  5,6  6,0  6,4  km/s 


Fig.  8.10. 


P  wave  velocity  versus  depth  for  a  medium  granite  after  BIRCH  (1958) 
according  to  a  temperature  gradient  of  8  °C/km. 


8.2.6  Depth  of  the  Mohorovicic  discontinuity 

The  Mohorovicic  discontinuity  (short:  Moho)  defines  the  crust-mantle  boundary;  it  is 
characteristically  a  rather  abrupt  increase  of  vp  from  less  than  7  km/s  to  values  around  8.0 
km/s  (extreme  values  of  7.5  and  8.5  km/s  are  rare).  The  :ncic''se  is  not  truly  a  first  order 
discontinuity  although  for  long  wavelengths  it  may  be  modelled  that  way;  in  reality  it  will 
always  be  a  transition  zone  between  crust  and  mantle  of  up  to  several  kilometers  vertical 
extent.  Therefore  it  is  a  matter  of  definition  at  which  depth  exactly  the  Moho  lies  and  one 
must  pay  attention  to  this  if  one  wishes  to  accept  published  results.  A  number  of  possible 
definitions  has  been  discussed  by  MEISSNER  (1986).  According  to  him  the  most  common 
definition  is  the  top  of  the  transition  zone  (or  layer)  of  the  greatest  vertical  P  wave  velocity 
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Fig.  8.12.  Tectonic  sketch  map  for  the  Swiss  Alps  taken  from  KL1 
(1982).  (1)  Molasse  sediments,  (2)  +  (3)  Swiss  jura,  (4)  -i 
Penninicum,  (7)  Eastern  Alpin,  (8)  Southern  Alpin ,  (9)  E 
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gradient  dvp/dz  leading  to  values  of  7.5  to  8.5  km/s.  This  definition  is  not  meant  to  apply 
to  other  high-velocity  gradient  zones  (MOSTAANPOUR  1984)-  Some  authors  define  the 
Moho  depth  at  the  bottom  of  the  transition  zone  or  at  its  center.  It  is  also  quite  common  to 
define  the  Moho  depth  Zra  by  the  exact  velocity  value  of  8.0  km/s.  Keeping  these 
uncertainties  of  definition  in  mind,  the  uncertainty  of  the  depth  Zra  for  the  Alps  has  been 
estimated  by  various  authors  to  be  ±  2  km  (EGLOFF  1979),  ±  1  ...  2  km  (PERRIER, 
RUEGG  1978)  and  ±  5  km  (MOSTAANPOUR  1984). 

Since  the  Moho  attracts  greatest  interest  there  are  more  data  on  it  than  on  any  other 
boundary  and  the  spatial  resolution  is  greatest.  In  this  work  a  number  of  more  detailled 
maps  and  crustal  sections  have  been  combined  with  more  regional  Moho  models  (e.g., 
MOSTAANPOUR  1984;  MEISSNER  1987).  Because  of  the  numerous  data  points  a  routine 
was  used  where  the  interpolating  function  is  a  fifth-degree  polynomial  within  each  triangle 
of  the  x-y  surface  of  the  data.  In  order  to  avoid  edge  effects  of  this  method  it  was  necessary 
to  include  a  0.75°  to  1°  strip  outside  the  actual  study  area.  The  result  is  shown  in  Fig.  8.13. 

Multiple  coverage  of  some  regions  from  several  publication  also  allows  us  to  estimate  the 
confidence  intervals  of  Zm  in  an  empirical  manner.  The  results  are  not  given  as  usual  as  a 
symmetric  standard  deviation  (±  AZ)  but  asymmetrically  for  the  shallowest  and  the 
deepest  Moho  depth  for  each  6’  *  10’  quadrangle.  For  the  best  cases  the  minimum  uncer¬ 
tainty  was  assumed  to  be  ±  1  km.  In  some  regions,  however,  we  found  differences  of  Zm  > 
10  km.  If  those  discrepancies  could  be  traced  to  the  differences  of  Moho  definition  in  the 
various  publications,  the  data  were  excluded  from  the  interpolation,  but  they  were  still  used 
to  estimate  the  confidence  limits.  Discrepancies  caused  by  unknown  effects  were  averaged  or 
treated  as  described  above. 

For  the  arc  of  the  western  Alps  and  for  the  Swiss  Alps  there  were  about  100  discrete  data 
points  available  ( CLOSS ,  LABROUSTE  1968;  KISSLING  1982).  In  addition  there  were  a 
few  published  contour  maps  for  these  regions.  The  data  had  to  be  fitted  together  mainly  in 
southwestern  Switzerland.  Here  also  contour  maps  of  the  whole  Alpine  region  (GIESE  et  al 
1976)  were  taken  into  account.  The  Alpine  forelands  in  the  north  and  south  were  represent¬ 
ed  only  by  contour  maps  (DRISLER,  JACOBY  1988;  GIESE  et  al  1976;  BARAZANGI, 
BROWN  1986)  which  are  consistent  with  other  representations.  Only  in  the  regions  of  the 
Swiss  Jura  and  the  southern  Rhine  Graben  we  noticed  significant  discrepancies  in  the  Moho 
depth  The  eastern  part  of  Austria  is  not  well  known  with  only  30  data  points  along  three 
seismic  profiles  (ARIC  1987)  available.  The  Jugoslavian  part  of  the  study  area  is  not 
covered  by  reliable  data;  we  had  to  revert  to  extrapolation.  The  data  density  is  much 
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greater  in  the  western  part  of  the  test  area  than  in  the  eastern  part,  and  the  data  sources 
are  also  more  detailled  and  reliable  in  the  west. 

The  upper  and  lower  confidence  limits  of  Zm  are  represented  as  two  sets  of  discrete  integer 
values  of  length  in  kilometer.  They  are  shown  in  Figs.  8.14  and  8.15.  A  representation  in 
the  form  of  contour  lines  is  of  no  meaning  and  no  value. 

The  contour  map  of  the  Mohorovicic  discontinuity  clearly  shows  the  arc  structure  of  the 
Alps  The  Zra  values  continuously  increase  from  about  32-38  km  in  the  region  of  the 
northern  Molasse  across  the  strike  direction  towards  the  axis  of  the  Alps  to  40-45  km  in  the 
Helveticum  and  Peninnicum.  In  West  Alpine  maximum  Zm  values  of  nearly  60  km  are 
reached.  In  contrast,  the  crust-mantle  boundary  in  the  central  East  Alpine  is  only  45-50 
km  deep.  In  the  Po  basin  the  crustal  thickness  is  only  25-35  km,  distinctly  less  than  in  the 
Molasse  to  the  north.  This  is  similar  to  the  sedimentary  thickness. 

Our  contour  map  (Fig.  8.13)  differs  clearly  from  other  interpretations  (e.g.,  MOSTAAN- 
POUR  1984;  MEISSNER  1987)  in  that  way  that  Zm  gradually  varies  from  the  center  of  the 
Alps  to  the  south.  The  maps  mentioned  show  a  rise  of  the  crust-mantle  boundary  from  the 
Po  basin  to  the  margin  of  the  western  Alpine  arc  to  about  10  km  from  which  the  boundary 
abruptly  drops  to  values  of  about  50  km  to  the  north  of  the  Insubric  Line.  This  behavior  is 
an  expression  of  the  Ivrea  body,  which  is  also  evident  in  the  partial  overlap  of  the  contour 
lines  manually  constructed.  Such  a  small-scale  structure  cannot  be  resolved  with  the  chosen 
grid,  and  the  interpolation  procedures  can  only  produce  smooth  surfaces  with  no  disconti¬ 
nuities.  The  Ivrea  zone  has  the  effect  that  the  contours  are  quite  dense  in  this  zone. 

The  Z,„  values  of  the  Po  plain  belong  to  the  "Adriatic  plate".  Crustal  doubling  as  postula¬ 
ted,  e  g.,  by  WIGGER  (1984)  and  GIESE  (1985)  cannot  be  recognized  in  our  model  and 
was  not  taken  into  account  in  our  estimate  of  the  confidence  limits.  There  is  no  evidence  for 
the  Eurasian  plate  to  have  been  pushed  a  long  distance  south  of  the  Insubric  line  below  the 
Adriatic  plate.  Locally  at  the  contact  of  the  Eurasian  and  Adriatic  plates  there  is  seismic 
evidence  for  crustal  subduction  at  the  southern  margin  of  the  west  Alpine  arc  and  the 
coastal  regions  of  Toscana  till  the  island  Elba.  A  continuous  connection  between  these  two 
regions  could,  however,  not  be  verified  by  refraction  seismology. 

Finally  we  note  the  shallow  Moho  with  Z„,  values  of  22  ...  25  km  below  the  southern  Rhine 
Graben.  However,  in  contrast  to  seismic  basement  there  is  no  abrupt  change  towards  the 
Black  Forest  and  Vosges  massifs. 
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Fig.  8.14.  Uncertainty  of  the  depth  of  the  Mohorovicic  discontinuity  in  {km]  in  case  of  a 
possible  lower  crust/mantle  boundary. 
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Fig.  8.15;  Uncertainty  of  the  depth  of  the  Mohorovicic  discontinuity  in  [km]  in  case  of  a 
possible  deeper  crust/mantle  boundary. 
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The  possible  deviations  of  the  Moho  depth  from  the  values  in  Fig.  8.13  are  in  large  areas 
less  than  ±  2  km  (Figs.  8.14,  8.15).  Possibly  greater  Zm  values  may  be  suspected  in  two  re¬ 
gions:  in  the  Ivrea  zone  to  the  southern  end  of  the  west  Alpine  arc  and  broadly  in  the 
eastern  Alps.  The  main  reasons  are  the  problems  of  incorporating  a  complicated  structure 
as  that  of  the  Ivrea  zone  and  considerable  discrepancies  in  different  Moho  maps  for  Europe 
(MEISSNER  1987)\  here  the  50  km  contour  extends  about  100  km  further  east  than  in  all 
other  maps.  Regions  with  possibly  smaller  Zra  values  are  again  along  the  west  Alpine  arc 
and  especially  in  the  Ivrea  Zone,  as  well  as  in  the  central  East  Alpine.  Also  the  Swiss  part 
of  the  northern  Molasse  basin  shows  uncertainties  towards  shallower  depths  by  3  ...  5  km. 
In  Figs.  8.14  and  8.15  it  is  evident,  that  generally  (with  the  exception  of  the  Ivrea  zone)  the 
reliability  and  accuracy  in  the  western  (Swiss)  part  of  the  study  area  are  much  better  than 
in  the  east.  The  reasons  are  the  higher  data  density  and  the  better  seismic  model. 


8.2.7  Thickness  of  the  lithosphere 

The  lithosphere  plays  an  important  geodynamic  role  as  the  uppermost  quasi-rigid  moving 
plate.  It  consists  of  the  crust  and  the  uppermost  part  of  the  mantle.  The  transition  to  the 
subjacent  asthenosphere  layer  is,  though  gradual,  a  drastic  change  in  rheology,  i.e.,  a 
decrease  of  viscosity  by  several  orders  of  magnitude.  The  asthenosphere  is  also  a  layer  of 
low  seismic  velocities  which  may  be  explained  by  partial  melt.  However,  the  transition  is 
principally  not  resolvable  as  a  reflecting  or  refracting  boundary.  The  depth  Zi  can  only  be 
determined  with  the  aid  of  surface  wave  dispersion  analysis  and  by  teleseismic  travel  time 
residuals.  Both  will  never  allow  the  definition  of  a  sharp  boundary.  Another  indirect 
method  to  determine  the  thickness  of  the  lithosphere  is  by  modelling  the  terrestial  heat  flow 
(tERMkK  1982).  However,  this  kind  of  modelling  is  ambiguous  and  depends  on  many 
vague  assumptions. 

The  present  map  (Fig.  8.16)  was  based  on  published  work  by  BABU§KA  et  al  (1985), 
SUHADOLC,  PANZA  (1988),  and  PANZA  et  al  (1980).  PANZA  (1981)  has  estimated  the 
uncertainty  to  be  ±  30  km.  Fig.  8.16  differs  in  parts  considerably  from  other  published 
maps.  For  example,  PANZA,  MULLER  (1980)  postulated  a  lithospheric  root  below  the 
Alps  up  to  200  km  thickness  and  thermally  induced  high  seismic  velocities  on  the  basis  of 
travel  time  residuals  (particularly  in  the  east  and  below  the  west  Alpine  arc).  This  root  has 
not  been  included  by  GEISS  (1987),  but  he  assumes  a  generally  thickened  lithoshere  (to 
130  km)  below  the  Alps.  At  present,  the  existence  of  the  lithosphere  is  being  debated  alto¬ 
gether,  e.g,  by  VIEL  (1987)  on  the  basis  of  observations  of  travel  time  residuals  at  his  Al- 


pine  stations  and  of  model  computations.  The  residuals  were  explained  by  him  by  crustal 
thickening. 

Fig.  8.16  represents  the  version  most  acceptable  to  us.  There  are  two  thickness  maxima 
near  the  west  Alpine  arc  and  below  the  central  eastern  Alps.  Towards  the  margins  of  the 
Alps  the  lithospheric  thickness  decreases  strongly  to  values  of  110  ...  130  km.  A  special  fea¬ 
ture  is  the  southern  Rhine  Graben  with  a  thickness  of  only  90  ...  100  km.  But  here  the  un¬ 
certainty  probably  greatly  exceeds  the  value  of  ±  30  km,  particularly  in  view  of  the  45  ... 
50  km  given  by  PANZA  et  al  (1980)  for  the  whole  European  rift  system.  The  same  is  true 
for  the  transition  from  the  Alps  to  the  Dinarides  in  the  southeastern  part  of  the  study  area, 
however,  with  the  opposite  sign.  PANZA  et  al  (1980)  give  values  >  120  km  for  that  area. 
Generally  one  has  to  take  into  account  that  surface  wave  dispersion  analysis  cannot  give 
much  spatial  resolution  horizontally  since  dispersion  can  only  become  evident  along  a  cer¬ 
tain  horizontal  path  length;  therefore  features  of  small  lateral  extend  will  appear  smoothed 
if  recognizable  at  all. 


8.2.8  Perturbations  of  seismic  velocities 

The  interpretation  is  based  on  the  same  data  as  the  construction  of  the  contour  map  of  Z|> 
does  Digitizing  the  contour  maps  was  difficult,  since  the  manually  constructed  contour 
lines  of  MOSTAANPOUR  (1984)  are  incomplete  and  appear  topologically  incorrect  at 
places  Therefore  we  took  more  of  the  original  seismic  profiles  into  account.  In  a  broad 
sense,  however,  the  results  are  comparable  to  those  of  MOSTAANPOUR  (1984).  We  deter¬ 
mined  the  lateral  variation  of  the  mean  P  wave  velocity  in  certain  regions  of  the  crust,  i.e., 
from  the  depth  range  from  the  surface  to  the  seismic  basement:  vp  (Fig.  8.17),  and  for  the 
whole  crust:  vm  (Fig.  8.18).  From  these  values  one  can  compute  the  mean  velocity  Vic  of 
the  crust  below  the  sediments,  i.e.,  from  Z\,  to  Zm  (Fig.  8.19);  according  to  MOSTAAN¬ 
POUR  (1984) 

AZ 

vie  =  with  AZ  =  Zm  -  Zb  .  (8— 3) 

vm  Vb 

Finally  the  Pn  velocity  vcm  at  the  crust  mantle  boundary  has  been  digitized.  It  is  the 
velocity  immediately  below  the  Mohorovicic  discontinuity  (Fig.  8.20). 
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Thickness  of  lithosphere  derived  from  seismological  observations  in  the  Alps. 
Contour  interval  10  km. 


10 


Fig. 


.17.  Mean  seismic  P  wave  velocity  in  the  upper  crust  vb  (layer  from  surface  down 
to  seismic  basement).  .Contour  interval  0.1  km/s. 
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apn?  n»! 


Seismic  velocity  vcra  at  the  crust/mantle  boundary  (uppermost  part  of  the 
mantle).  Contour  interval  0.1  km/s. 


The  contours  of  vb  again  show  the  typical  structure  of  the  Alps.  The  Molasse  basins  and  the 
southern  Rhine  Graben  filled  with  sediments  are  characterized  by  low  seismic  velocity  as 
theoretically  expected.  In  the  central  East  Alpine  and  the  western  Alpine  arc  vb  reaches  re¬ 
latively  high  values  of  5.6  ...  5.8  km/s.  It  may  be  suspected  that  the  definition  of  the  seis¬ 
mic  basement  (see  section  8.2.5)  ist  not  valid  there.  The  variation  of  the  mean  P  wave  ve¬ 
locity  vm  for  the  whole  crust  is  much  smaller  (about  5.8  ...  6.3  km/s)  than  that  within  the 
upper  crust.  The  variation  correlates  well  with  the  locations  of  the  sedimentary  basins.  We 
note,  however,  the  anomalously  low  average  velocity  in  the  western  Alps  (5.8  km/s).  Here 
the  mean  P  wave  velocity  is  the  same  in  the  upper  and  lower  crustal  layers.  This  is  also 
noticeable  in  the  mean  velocities  vic  in  the  sediment-stripped  crust  (Fig.  8.19).  Otherwise 
the  mean  velocities  in  the  lower  crustal  regions  range  between  6.0  and  6.5  km/s.  The 
maximum  values  are  found  below  the  sedimentary  basins.  This  agrees  with  results  of 
BIRCH  (1958). 

Note  that  Fig.  8.19  clearly  differs  from  the  picture  presented  by  MOSTAANPOUR  (1984). 
The  reason  for  this  is  partly  that  we  computed  the  average  velocity  vic  for  each  6'  x  10' 
quadrangle,  while  MOSTAANPOUR  (1984)  had  computed  v i,-  only  for  a  few  points  and 
constructed  the  contour  lines  manually. 

The  Pn  velocity  vrm  just  below  the  crust-mantle  boundary  shows  relatively  low  values 
(7.8  km/s)  in  the  region  of  the  deepest  Moho.  These  regions  are  surrounded  by  zones  of  re¬ 
latively  high  Pn  velocity  (8.2  km/s).  There  is  some  north-south  asymmetry  in  the  sedimen¬ 
tary  basins;  the  Molasse  is  underlain  by  high  velocities  (8.0  ...  8.2  km/s)  while  the  Po  plain 
is  underlain  by  relatively  low  velocities  (7.8  km/s).  The  anomalously  low  velocities  in  the 
lower  crust  of  the  western  Alps  has  no  counterpart  below  the  crust  mantle  boundary,  but 
there  are  anomalously  low  Pn  velocities  below  the  southern  Rhine  Graben  (7.5  ...  7.6  km/s). 


8.2.9  Digital  density  model 

In  order  to  compute  topographic  (terrain)  reduction  of  gravity  measurements,  one  should 
know  the  density  distribution  near  the  earth’s  surface  above  the  reference  level.  For  compu¬ 
ter-adequate  procedures  one  needs  a  digital  density  model  besides  the  digital  elevations. 
The  density  surface  rocks  range  between  2.0  (unconsolidated  sediments)  and  3.0  g/cm3  (ul- 
trabasic  magmatic  n><  k s )  The  usual  assumption  of  an  average  density  of  2  67  g/cm3  for  the 
;,ui  be  ml  <  r  i  tal  io'  k  .  I<;ep.  to  un"ii,iml  i<s  or  i-ii'irs  m  He  <|e|<  i  mi  nation  ol  tie-  topogiaph 
ic  reductions.  For  Austria  a  digital  density  model  is  now  available  (see  Fig.  8.21).  It  grew 
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Distribution  of  Surface  Density 


out  of  combination  of  data  sets  already  published  and  a  special  data  collection  for  the  pres¬ 
ent  purpose.  The  bases  are  well  researched  densities  of  surface  rocks  in  Austria.  With  the 
aid  of  a  geological  map  the  region  was  divided  into  density  provinces.  This  regionalisation 
together  with  the  corresponding  density  values  is  shown  in  Fig.  8.22.  The  densities  are 
given  to  the  nearest  0.05  g/cm3  and  range  from  2.0  to  2.85  g/cm3.  WALACH  (1987)  gives 
the  standard  error  of  his  model  as  4  %  and  assumes  that  this  holds  from  the  earth’s  surface 
to  the  sea  level.  These  values  had  essentially  been  determined  from  samples  by  weighting  in 
air  and  water,  by  applying  Nettleton’s  gravimetric  method,  and  by  gravity  measurements 
in  vertical  mine  shafts.  » 


Fig.  8.22.  Map  of  surface  rock  densities  in  [g/cm3]  in  Austria  according  to  STEINHAU- 
SER  et  al  (l98S). 
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9. 


NUMERICAL  INVESTIGATIONS 
9.1  Prediction  of  gravity  anomalies  in  test  area  "Rossdorf" 

The  purpose  of  the  computations  was  to  investigate  to  what  extent  (surface)  density  data 
can  improve  the  prediction  accuracy  of  gravity  anomalies.  This  includes  also  the  case  where 
gaps  without  gravity  observations  can  be  overbridged  via  density  observations.  For  the  nu¬ 
merical  realization  methods  were  preferred  fitting  into  the  general  integrated  geodesy  ad¬ 
justment  scheme. 

From  the  various  models  discussed  in  chap.  4  two  were  implemented  in  FORTRAN77  pro¬ 
grams:  The  method  of  quasi-harmonic  inversion  and  the  attenuated  white  noise  gravity  cova¬ 
riance  model  of  JORDAN,  HELLER  (1978). 

For  the  investigations  the  stochastic  properties  of  the  gravity  anomalies  as  well  as  those  of 
the  density  data  had  to  be  evaluated.  For  the  fulfillment  of  the  condition  E{-}  =  0  mean 
values  as  simple  trend  were  subtracted.  Topographic-isostatic  reduction  of  gravity  anoma¬ 
lies  was  done  using  the  prism  integration  method  (FORSBERG  1984),  the  digital  terrain 
model  (100  *  100  m3)  and  the  compensation  model  of  Airy-Heiskanen  type  (HEISKANEN, 
MORITZ  1967,  p.  185),  however,  with  the  mean  surface  density  of  p  -  2.63  g/cm3  of  all  78 
stations.  For  the  determination  of  the  Bouguei  anomalies  the  corresponding  discrete  density 
data  of  the  stations  were  applied.  The  resulting  characteristics  of  the  covariance  functions 
are  summarized  in  Tab.  9.1. 

Table  9.1.  Characteristic  parameters  of  empirical  covariance  functions  in  test  area 
"Rossdorf" 


Parameters  of 
covariance  functions 

Free-air 

gravity 

anomalies 

AgF 

Bouguer 

gravity 

anomalies 

Ago 

Topographic- 
isostatic 
anomal ies 

AgTI 

Surface  1 

density 

data 

8p 

Variance 

144  mGal3 

93  mGal3 

103  mGal3 

0.014(g/cm3)3 

Correlation  length  ip 

5.9  km 

5.8  km 

5.7  km 

2.3  km 

Variance  of 
horizontal  gradient 

570  E3 

736  E3 

780  E3 

0.1(g/cm3)/km 
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As  a  second  step  an  iterative  fit  to  the  analytical  covariance  model  of  TSCHERNING, 
RAPP  (1974)  to  the  empirically  determined  covariance  functions  was  carried  out  (for  use  in 
the  quasi-harmonic  inversion  method).  The  results  can  be  found  in  Tab.  9.2. 


Table  9.2.  Parameters  of  the  TSCHERNING,  RAPP  (1974)  covariance  model  for  the 
test  area  "Rossdorf"  using  Agji. 


Order  of  local  covariance  function  n 

619 

(Rb/Re)2 

0.999615 

AVAR 

96.91 

In  Tab.  9.2  (Rd/Re)2  is  the  square  of  the  ratio  of  the  radius  Rn  of  the  so-called  Bjerham- 
mar  sphere  to  the  earth’s  radius  Re,  and  AVAR  is  a  scale  factor  of  the  model.  Comparison 
of  the  model  crosscovariances  cov(Ag ,6p)  with  the  corresponding  empirical  values,  however, 
showed  a  significant  disagreement  (Figs.  9.1  to  9.3) 

The  third  step  was  the  prediction  of  gravity  anomalies.  The  method  applied  was  the  follow¬ 
ing:  We  assumed  for  each  point  Pi  no  observation  being  available  and  predicted  for  Pi  the 
gravity  anomaly  from  its  77  (=  78  -  1)  neighbours.  The  difference  between  the  available 
"observation"  Ag  and  the  predicted  value  Agprej  can  be  considered  as  the  contribution  of 
point  Pi  to  a  "true"  prediction  error  defined  by 


a  = 


t 


i  =  i 


(Ag  —  Agpred)i  /  n. 


(9-1) 


The  predictions  were  done  using  gravity  anomalies  Ag  alone,  a  combination  of  gravity  ano¬ 
malies  and  surface  densities,  and  also  using  density  values  alone. 


Prediction  results  using  the  quasi-harmonic  inversion  method 
Results  of  the  test  runs  are  given  in  Tab.  9.3. 
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Table  9.3.  Results  of  gravity  anomaly  prediction  in  the  test  area  "Rossdorf1  using  the 
quasi-harmonic  inversion  method.  "True"  prediction  errors  in  [mGal]. 


Data 

AgF 

AgB 

agTi 

Gravity  anomalies  alone 

2.9 

2.8 

Gravity  anomalies  and 
density  data 

5.9 

4.4 

4.4 

Density  data  alone 

5.4 

4.3 

4.4 

The  prediction  of  gravity  anomalies  using  only  gravity  anomalies  was  more  or  less  carried 
out  in  order  to  get  some  kind  of  reference.  As  expected  the  consideration  of  a  simple 
Bouguer  plate  as  height  model  (-•  Age)  resulted  already  in  a  significant  improvement  in  the 
gravity  prediction.  Whereas  the  use  of  a  detailed  elevation  grid  assuming  in  addition  an 
isostatic  compensation  could  not  show  any  further  reduction  in  the  "true"  prediction  error. 
This  could  possibly  be  explained  that  for  the  determination  of  the  Bouguer  plate  effect 
observed  (point)  density  values  were  used  and  secondly,  that  it  might  be  questioned 
whether  in  that  area  the  consideration  of  an  isostatic  compensation  mechanism  is  justified. 
It  is  surprising  that  the  free-air  gravity  anomaly  prediction  from  density  data  alone  was 
only  slightly  worse  (±  5.4  mGal)  than  the  one  from  gravity  anomalies  (±  4.0  mGal). 

What  at  the  moment  cannot  reasonably  be  explained  is  the  fact  that  no  improvement  in 
gravity  prediction  was  achieved  when  using  a  combined  set  of  gravity  anomalies  and 
density  data. 

Gravity  prediction  using  the  attenuated  white  noise  gravity  covariance  model 

For  numerical  investigations  the  spherical  (6-19)  as  well  as  the  flat-earth  form  (6-21)  of  the 
attenuated  white  noise  gravity  covariance  model  (see  chap.  6.2)  were  programmed.  The  ex¬ 
pressions  for  the  necessary  covariance  and  crosscovariance  functions  can  be  found  in  Appen¬ 
dix  A. 

By  varying  the  parameters  D  (characteristic  depth)  and  Cot  (variance  of  the  disturbing  po¬ 
tential)  of  the  attenuated  white  noise  gravity  covariance  model  the  following  empirical 
relations  between  D,  Cot  and  ,  ip  (variance  and  correlation  length  of  the 

autocovariance  function  of  gravity  anomalies)  were  found: 
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(9-2) 


D  =  U, 


C 


OT 


TF 


1  + 


D 

K 


MT 


(9-3) 


or,  when  using  the  asymptotic  forms  (flat-earth  model), 


C«Ar  D2 


CoT  ~  15  DOC 


(9-4) 


The  used  units  are: 


D,  ip,  R  [km]  , 

Cot  [(m/s)4]  , 

C0Ag  fa^2! 


For  the  "Rossdorf  test  area  the  free  parameters  are  determined  as  function  of  the  variance 
and  the  correlation  length  of  the  autocovariance  function  of  the  corresponding  gravity 
anomalies  (free-air,  Bouguer,  topographic-isostatic). 


Table  9.4.  Characteristic  depth  D  and  variance  of  the  disturbing  potential  C0t  as 
function  of  the  correlation  length  ip  and  the  variance  C0Ag  of  the 
corresponding  gravity  anomaly  covariance  function. 


Free-air 

gravity 

anomalies 

AgF 

Bouguer 

gravity 

anomalies 

AgB 

Topographic- 

isostatic 

anomalies 

AgTI 

Variance  ^o^g 
Correlation  length  ip 

144  mGal2 

5.9  km 

93  mGal2 

5.8  km 

103  mGal2 

5.7  km 

Characteristic  depth  D 
Model  variance  Cot 

7.9  km 
0.60  (!)< 

7.8  km 
0.38  (|)4 

7.7  km 

0.40  (f)4 

Plots  of  the  empirical  as  well  as  of  the  model  covariance  functions  of  the  Tscherning/  Rapp 
and  the  attenuated  white  noise  gravity  covariance  model  of  Jordan/Heller  are  shown  for 
comparison  in  Figs.  9.1  to  9.3.  The  graphs  indicate  that  the  attenuated  white  noise 
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Autocovariance  function  of  free-air  gravity  anomalies 


Crosscovariance  function  between  densities  and  free-air  anomalies 


Autocovariance  function  of  surface  density  values 


—  Empirical 

-  - Tscherning/Rapp 
•  -  •  Jordan/Heller 


Fig.  9.1.  Covariance  functions  for  the  test  area  "Rossdorf"  (Free-air  gravity  anomaly 
fit). 
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Autocovariance  function  of  Bouguer  gravity  anomalies 


Distance  [km] 


Crosscovariance  function  between  densities  and  Bouguer  anomalies 


Oistance  [km] 
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Autocovariance  function  of  isostatic  gravity  anomalies 


Crosscovariance  function  between  densities  and  isostatic  anomalies 


Autocovariance  function  of  surface  density  values 


—  Empirical 

-  - Tscherning/Rapp 
•••  Jordan/Heller 


Fig.  9.3.  Covariance  functions  for  the  test  area  "Rossdorf"  (Topographic-isostatic 
gravity  anomaly  fit). 
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gravity  covariance  model  fits  better  to  the  empirical  autocovariance  functions 
^ SpSp  ’  a^ough  ‘ts  var'ance  C 0(6p,Sp)  is  still  too  high.  Both,  crosscovariances  C 
autocovariances  C^-  ^  are  decreasing  too  slowly  with  respect  to  small  distances. 


CA6Ag’ 

«/>Ag  Md 


The  numerical  prediction  trials  were  carried  out  in  the  same  way  as  outlined  before  with  re¬ 
gard  to  the  use  of  the  quasi-harmonic  inversion  method.  Only  observations  within  a  radius 
of  approx.  6  km  (corresponding  to  a  little  bit  more  than  the  correlation  length  of  the 
empirical  gravity  covariance  functions,  see  Tab.  9.1)  were  considered  in  the  pointwise 
gravity  prediction  algorithm.  The  results  are  given  in  Tab.  9.5. 


Table  9.5.  Results  of  gravity  anomaly  prediction  in  the  test  area  "Rossdorf"  using  the 
attenuated  white  noise  gravity  covariance  model  of  Jordan/Heller.  "True"  pre¬ 
diction  errors  in  [mGal], 


Data 

AgF 

Ago 

AgTI 

Gravity  anomalies  alone 

4.3 

2.9 

2.8 

Gravity  anomalies  and 
density  data 

- 

- 

- 

Density  data  alone 

fc= - 

5.8 

5.6 

5.4 

If  the  reader  compares  the  results  in  Tab.  9.5  with  the  ones  in  Tab.  9.3  no  principal 
differences  can  be  found  except  that  when  using  the  Jordan/Heller  covariance  model  in  the 
form  like  here  -  one-layer  crust  -  it  is  not  possible  to  consider  gravity  and  density  data  in  a 
combined  observation  set  due  to  ill-conditioning  and  singularities,  resp.,  in  the  flat  earth 
one-layer-crust  covariance  model  (see  the  comments  after  eq.  (6-28)). 

Nevertheless,  for  further  studies  in  gravity  prediction  using  density  data  the  application  of 
a  Jordan/Heller  multiple-layer  crust  model  considering  regional  isostatic  compensation  (for 
example  in  the  form  of  eqs.  (6-29)  to  (6-29d))  seems  to  be  promising,  since  it  contains  no 
singularities  or  ill-conditioning  and  can  be  widely  adapted  to  the  physical  reality  (in 
contrast  to  the  quasi-harmonic  method). 
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An  interesting  by-product  of  the  developed  gravity  prediction  algorithms  is  their  use  for  the 
determination  of  density  from  gravity  anomalies  (inverse  problem).  The  results  of  these 
computations  can  be  found  in  Tab.  9.6. 


Table  9.6.  True  errors  of  density  prediction  in  [g  cm"3] 
in  the  test  area  "Rossdorf” 


Data 

Quasi-harmonic 

inversion 

Tscherning/Rapp  model 

Jordan/Heller 
covariance  model 

Gravity  anomalies  alone 

0.10 

0.08 

Gravity  anomalies  and 

0.16 

density  data 

Density  data  alone 

0.07 

0.08  ! 

Although  the  results  here  are  not  representative,  the  accuracy  of  estimated  density  in  the 
range  of  ±  0.1  g  cm  3  is  remarkable. 


9.2  Prediction  of  gravity  anomalies  in  the  European  Alps 

Whereas  in  chapter  9.1  test  computations  were  carried  out  in  a  very  local  area,  the  inten¬ 
tion  here  was  to  repeat  those  ones  on  a  more  regional  scale.  For  that  reason  the  data  collec¬ 
ted  over  the  European  Alps  -  described  in  chap.  8.2  -  were  used. 

Empirical  covariance  functions  over  the  European  Alps  {6°  <  X  <  16° 5  ;  45°  <  <p  <  48°5) 

For  the  topographic-isostatic  gravity  anomalies  (see  chap.  8.2.2)  the  empirical  covariance 
function  was  determined  using  some  kind  of  homogeneous  thinning  procedure  for  the  data, 
similar  to  the  one  described  by  LANDAU  et  al  (1988).  The  residual  trend  was  removed 
using  a  simple  linear  regression  of  type 

Agn  =  a  +  b  H  (9-2) 
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where  H  denotes  the  height  and  a,b  are  the  parameters  which  were  determined  to  be 

a  =  46.7  mGal 
b  =  -  0.029  mGal/m  . 

For  the  determination  of  the  corresponding  empirical  surface  density  covariance  function  a 

simple  constant  (mean  value)  as  trend  was  considered  to  be  sufficient,  p  =  2.41  g/cm3.  The 
results  of  the  computations  are  summarized  in  Tab.  9.7. 


Table  9.7.  Characteristics  of  empirical  covariance  functions  over  the  European  Alps 
(6»  <  A  <  16 ?5  ;  45°  <  tp  <  48°5). 


Parameters  of 
covariance  functions 

Topographic-isostatic 
gravity  anomalies 

Agn 

Surface  density  data 
(3'  ■  5') 

6p 

Variance 

660  mGal2 

0.15  (g/cm3)2  | 

Correlation  length  ip 

24.3  km 

46.8  km 

Variance  of 

208  E2 

horizontal  gradient 

| 

Plots  of  the  empirical  covariance  functions  are  found  in  Figs.  9.4  and  9.5. 

Comparing  parameters  for  the  gravity  anomaly  covariance  function  with  the  results  of 
(SVNKEL  et  al  1983)  shows  a  difference  in  the  variance.  This,  however,  can  be  easily  ex¬ 
plained  by  the  fact  that  they  restricted  themselves  to  data  covering  Austria  only  whereas 
we  also  included  gravity  anomalies  for  the  rest  of  the  European  Alps.  Especially  the  gravity 
field  in  Switzerland  and  northern  Italy  shows  large  variations  in  the  gravity  anomalies, 
compare  Tab.  9.8. 


Table  9.8.  Variations  in  gravity  anomalies  in  [mGal]. 


Austria 

F.R.  Germany 

Italy 

Switzerland 

Minimum 

-  48 

-  60 

-  108 

-  60  | 

!  Maximum 

ii 

♦  121 

+  3 

+  160 

+  134 
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Covariances  [g/cm3]##2  ^  Covariances  [mgal]#*2 


Empirical  covariance  functions  over  the  1 0  *  1°  block  (13°  <  A  <  14°  and  47®  <  <p  <  483) 

In  order  to  save  computer  time  it  was  decided  to  test  the  quasi-harmonic  inversion  in  the 
1°  x  1°  block  mentioned  above.  There  the  topography  varies  from  500  m  up  to  2000  m  (Fig. 
8.4).  Within  this  area  an  amount  of  3396  gravity  values  was  found  in  the  data  base. 


Table  9.9.  Characteristics  of  empirical  covariance  functions  over  the  1°  *  l»  block 
(13°  <  A  <  14°  and  47°  <  ip  <  48°)  of  the  European  Alps. 


larameters  of 
covariance  functions 

Topographic-isostatic 
gravity  anomalies 

AgTI 

Surface  density  data 
(3'  x  5') 
bp 

Variance 

472  mGal2 

0.013  (g/cm3)2 

Correlation  length  ip 

18.5  km 

20.3  km 

Variance  of 

313  E2 

horizontal  gradient 

Similarly  to  the  computations  in  the  "Rossdorf"  test  area  an  iterative  fit  of  the  analytical 
covariance  model  of  TSCHERNING,  RAPP  (1914)  was  carried  out  which  resulted  in  the 
parameters  outlined  in  Tab.  9.10. 

Table  9.10  Parameters  of  the  TSCHERNING,  RAPP  (1974)  covariance  model  for  the 
1°  x  io  block  of  the  European  Alps. 


Order  of  local  covariance  function  n 

200 

(Rd/Re)2 

0.998615 

AVAR 

517.28 

Test  computations  in  gravity  interpolation  were  again  done  in  the  same  manner  as  in  the 
"Rossdorf"  area  using  the  quasi-harmonic  inversion  method.  The  results  are  found  in  Tab. 
9.11 
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Table  9.11  Results  of  gravity  anomaly  prediction  in  the  European  Alps  (1°  *  1°  block) 
using  the  quasi-harmonic  inversion  method. 


Data 

True  prediction  error 

Topographic-isostatic 

5.2  mGal 

gravity  anomalies 

Topographic-isostatic 

7.0  mGal 

gravity  anomalies 

and  density  data 

Density  data  alone 

failed 

The  results  have  to  be  interpreted  as  "failed  with  respect  to  the  assumption  of  quasi- 
harmonicity  for  density".  This  is  also  in  agreement  with  results  published  by  HEIN  et  al 
(1988)  and  experiences  of  STRYKOWSKI,  TSCHERNING  (1989,  pers.  comm.)  who  also 
got  from  numerical  investigations  the  impression  that  the  "quasi-harmonicity"  is  a  too 
strong  condition  when  using  it  to  solve  the  inverse  gravitational  problem.  There  might  be, 
however,  cases  where  the  method  gives  reasonable  results  with  respect  to  gravity  prediction. 
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10.  CONCLUSIONS 


The  study  presents  some  first  trials  to  use  density  and  seismic  data  in  a  direct  way  in  grav¬ 
ity  field  approximation,  in  particular  gravity  prediction.  Although  the  geodesists  are 
primarily  not  interested  in  the  inverse  problem,  it  is  necessary  to  deal  with  it  to  some 
extent.  The  geophysicist’s  view  on  the  structure  of  the  earth’s  crust  and  its  knowledge 
through  seismics  can  help  to  improve  the  geodetic  gravity  prediction  algorithms.  In 
particular,  empirical  relationship  between  seismic  velocity  and  density  might  be  a  valuable 
replacement  of  complex  wave  propagation  formulas  (see  chapter  3). 

In  the  theoretical  part  a  few  new  approaches  are  presented  which  seem  to  be  very  promis¬ 
ing.  Since  methods  developed  up  till  now  use  rather  mathematical  than  (realistic)  physical 
constraints  for  the  covariance  propagation  needed  in  prediction  algorithms  of  collocation 
type,  special  emphasis  was  put  to  take  advantage  of  Newton’s  attraction  integral  for  that 
purpose.  The  generalized  isostatic  response  theory  of  DORMAN,  LEWIS  (1970)  seems  to  be 
a  promising  candidate  in  the  gravity  prediction  applications. 

Although  the  quasi-harmonic  inversion  is  based  on  physically  not  realistic  constraints,  it 
may  offer  in  some  cases  reasonable  prediction  results.  HELLER’s  1976  developed 
attenuated  white  noise  statistical  gravity  model  seems  to  be  more  realistic  when 
introducing  multiple  layers  of  the  earth’s  crust.  In  particular,  the  proposal  of  JORDAN 
(1978)  to  extend  this  covariance  theory  to  regional  isostatic  compensation  might  result  in 
an  even  more  realistic  model  for  the  terrain  and  crust-mantle  contrasts. 

In  the  consideration  of  seismic  data  for  gravity  prediction  (chapter  7)  the  work  of  the 
geophysicist  CHERNOV  (1960)  is  appreciated.  His  statements  about  density  and  the  elastic 
Lame  parameters  at  being  randomly  distributed  in  space  forming  a  stochastic  process 
allowed  the  incorporation  of  density  and  seismic  data  in  the  integrated  geodesy  adjustment 
model  and  the  formulation  of  covariance  functions  according  to  seismic  wave  motion  in  an 
inhomogeneous  medium. 

The  detailed  numerical  investigations  in  chapter  8  in  a  local  as  well  as  in  the  European  alps 
brought  more  insight  in  the  statistical  behaviour  of  density  and  present  various 
characteristics  of  covariance  functions.  Some  first  gravity  predictions  using  density  and  the 
quasi -harmonic  inversion  method  as  well  as  the  attenuated  white  noise  statistical  gravity 
covariance  model  show  reasonable  results. 
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Since  this  report  can  represent  only  some  initial  work  in  this  field,  it  is  recommended  to 
continue  this  study,  in  particular,  by  realizing  numerically  (i)  the  proposed  covariance 
model  in  chapter  6.3  in  connection  with  the  generalized  isostatic  response  theory,  (ii)  the 
extension  of  HELLER’s  covariance  model  to  regional  compensation  as  proposed  by 
JORDAN  (1978),  and  (iii)  the  inclusion  of  seismic  data  based  on  the  statistical  assumptions 
found  by  CHERNOV  (1960). 

There  is  no  doubt  that  density  and  seismic  information  will  strengthen  a^nd  improve  gravity 
prediction,  in  particular  in  areas  where  gravity  gaps  are.  However,  the  best  method  still  has 
to  be  found  in  the  future.  It  is  hoped  that  this  report  provides  a  platform  for  further  work. 


APPENDIX  A 


Auto-  and  crosscovariances  of  the  gravity  field  functionals  based  on  the  attenuated  white  noise 
statistical  gravity  model 


By  HELLER  and  JORDAN  (1979)  the  autocovariance  function  of  the  disturbing 
potential  T  is  given  in  the  form 

2  2  2  2  4 

D  (2R-D)  [  rl  r2  -(R-D)  ]  Co(T) 

Cov(T.T) (rl ,r2 ,Psi)  -  . - . 

2  2  4  2  2  2  3/2 

(2R  -2DR+D  )  [  (R-D)  +  rl  r2  -  2rlr2(R-D)  cos(Psi)  ] 


Substituting  in  (6-19)  D/2  by  D  and  introducing  the  abbreviations 

3  3  4  4 

A  -  D  ( 2R-D)  /  [  R  -  (R-D)  ] 

2  2  2  2 

-  D  (2R-D)  /  [  R  +  (R-D)  ] 

2  2  2  2  2 

-  [  R  -  (R-D)  ]  /  [  R  +  (R-D)  J 


2  2  4 

B  -  rl  r2  -  (R-D) 


2  2  4  2 

C  -  rl  r2  +  (R-D)  -  2rlr2(R-D)  cos(Psi) 


leads  to 

B 

Cov(T ,T) (rl , r2 , Psi)  -  A  .  Co(T)  . 

3/2 

C 

The  auto-  and  crosscovariances  between  the  functionals  of  the  gravity  field  were 
computed  with  the  help  of  the  Algebraic  Computation  Software  REDUCE  (RAYNA, 1987) . 
So  the  derivatives  of  T  and  Cov(T,T)  with  respect  to  rl,  r2,  ...  are  written  in 
the  form 

DF(T.rl) ,  DF(T,r2) ,  ... 

DF(Cov(T,T) , rl) . 

To  get  higher  derivatives  you  have  to  add  further  arguments,  e.g. 
DF(Cov(T,T),rl,2,r2) 

means  the  second  derivative  with  respect  to  rl,  first  with  respect  to  rl. 
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The  functionals  of  the  gravity  field  in  F  are  introduced  by  the  following 
expressions: 


1)  geoldal  height  N1 

2)  gravity  disturbance  d  gl 

3)  gravity  anomaly  D  gl 

4)  radial  gravity  gradient 

5)  second  radial  derivative  of  T 

€)  N-S -component  of  the 

deflection  of  the  vertical  Xil 

7)  E-W- component  of  the 

deflection  of  the  vertical  Etal 

8)  density  contrast  d  rho 


-  Tl/G 

-  -  DF(Tl,rl) 

-  -  DF(Tl,rl)  -  2*Tl/rl 

-  -  DF(Tl,rl,2)  -  2*DF(T1 , rl)/rl 

-  +  DF(Tl,rl,2) 

-  -  DF(T1 ,pl)/ (G*R) 

-  -  DF(T1 , ll)/(G*R*cos(pl) ) 

-  -  DF(Tl,rl)*DIRAC(rl  -  R’)/(4*Pi*k) 


In  the  formulas  above  G  is  the  normal  gravity  value,  k  the  Newtonian  gravity 
constant,  R  the  earth's  mean  radius,  R'  the  radius  of  the  density  contrast, 
pi  the  geographical  latitude  and  11  the  geographical  longitude  of  P. 

To  get  the  functionals  in  Q  substitute  Tl, rl.pl, 11  by  T2,r2,p2,12. 

The  covariances  between  the  functionals  are  computed  as  matrix  elements 
Cov(i.j),  where  i  and  J  are  the  numbers  of  the  functionals  in  P  and  Q, 
e.g.  Cov(3,8)  is  the  crosscovariance  between  the  gravity  anomaly  in  P  and 
the  density  contrast  in  Q. 


2 

Cov(l , 1)  -  Cov(T,T)/G 

Cov(l , 2)  -  -  DF(Cov(T,T) ,r2)/G 

Cov(l,3)  -  (  -  DF(Cov(T,T) ,r2)  -  2*Cov(T,T)/r2)/G 

Cov(l ,4)  -  (  -  DF(Cov(T,T),r2,2)  -  2*DF(Cov(T,T) ,r2)/r2)/G 

Cov(l , 5)  -  DF(Cov(T,T),r2,2)/G 

2 

Cov(l,6)  -  -  DF(Cov(T,T) ,p2)/(G  *R) 

2 

Cov(l , 7)  -  -  DF(Cov(T,T) , 12)/(G  *R*cos(p2)) 


Cov(l , 8)  -  -  DF(Cov(T,T) , r2)*DIRAC(r2  -  R'),  (4*Pi*k*G) 


Cov(2 , 1)  -  -  DF(Cov(T,T) ,rl)/G 
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Cov( 2,2) 


DF(Cov(T,T) ,rl,r2) 


Cov(2 , 3) 
Cov( 2,4) 
Cov(2 , 5 ) 
Cov( 2,6) 
Cov( 2 , 7 ) 
Cov(2 ,8) 

Cov( 3,1) 
Cov( 3,2) 
Cov( 3,3) 

Cov( 3,4) 

Cov( 3 . 5) 
Cov(3 ,6) 
Cov( 3,7) 
Cov( 3 , 8) 

Cov(4 , 1) 
Cov(4 , 2) 
Cov(4 , 3) 

Cov(4,4) 

Cov(4 , 5) 
Cov(4 , 6) 
Cov(4 , 7) 
Cov(4 , 8) 


DF(Cov(T,T) ,rl,r2)  +  2*DF(Cov(T,T) ,rl)/r2 
DF(Cov(T , T) , rl , r2 , 2)  +  2*DF(Cov(T,T) ,rl,r2)/r2 

-  DF(Cov(T , T) , rl , r2 , 2) 

DF(Cov(T,T) , rl ,p2)/(G*R) 

Dr(Cov(T , T) , rl , 12)/(G*R*cos (p2) ) 

DF(Cov(T,T) , rl , r2)*DIRAC(r2  -  R’)/(4*Pi*k) 

(  -  DF(Cov(T,T) , rl)  -  2*Cov(T,T)/rl)/G 
DF(Cov(T,T) ,rl,r2)  +  2*DF(Cov(T ,T) , r2)/rl 
DF(Cov(T,T) ,rl,r2)  +  2*DF(Cov(T,T) , rl)/r2 
+  2*DF(Cov(T,T) , r2)/rl  +  4*Cov(T,T)/(rl*r2) 

DF(Cov(T , T) , rl , r2 , 2)  +  2*DF(Cov(T,T) ,rl,r2)/r2 
+  2*DF(Cov(T,T) ,r2,2)/rl  +  4*DF(Cov(T,T) ,r2)/(rl*r2) 

-  DF(Cov(T , T) , rl , r2 , 2)  -  2*DF(Cov(T,T) , r2 , 2)/rl 
(DF(Cov(T,T) ,rl,p2)  +  2*DF(Cov(T,T) ,p2)/rl)/(G*R) 

(DF(Cov(T , T) , rl , 12)  +  2*DF(Cov(T,T) , 12)/rl)/(G*R*cos(p2) ) 

DIRAC(r2  -  R')*(DF(Cov(T,T),rl,r2)  +  2*DF(Cov(T,T) ,r2)/rl)/(4*Pi*k) 

(  -  DF(Cov(T,T) ,rl,2)  -  2*DF(Cov(T,T) , rl)/rl)/G 
DF(Cov(T,T) ,rl ,2,r2)  +  2*DF(Cov(T,T) , rl , r2)/rl 
DF(Cov(T,T) , rl,2,r2)  +  2*DF(Cov(T,T) ,rl,r2)/rl 

+  2*DF(Cov(T ,T) , rl , 2)/r2  +  4*DF(Cov(T,T) , rl)/(rl*r2) 

DF(Cov(T ,T) , rl , 2 , r2 , 2)  +  2*DF(Cov(T,T) ,rl,r2,2)/rl 
+  2*DF(Cov(T,T) ,rl,2,r2)/r2  +  4*DF(Cov(T,T) , rl , r2)/(rl*r2) 

-  DF(Cov(T ,T) , rl , 2 , r2 , 2)  -  2*DF(Cov(T,T) ,rl,r2,2)/rl 

(DF(Cov(T,T) , rl , 2 ,p2)  +  2*DF(Cov(T,T) , rl ,p2)/rl)/(G*R) 

(DF(Cov(T,T) , rl , 2 , 12)  +  2*DF(Cov(T,T) , rl , 12)/rl)/(G*R*cos(p2) ) 

(DIRAC(r2  -  R' )*(DF(Cov(T,T) , rl,2,r2)  +  2*DF(Cov(T,T) ,rl ,r2)/rl)/ 
(4*Pi*k) 
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Cov(5 , 1) 
Cov(5 , 2) 
Cov(5 , 3) 
Cov(5 ,4) 
Cov(5 , 5) 
Cov(5,6) 
Cov(5,7) 
Cov(5,8) 


-  +  DF(Cov(T,T) ,rl,2)/G 

-  -  DF(Cov(T,T),rl,2,r2) 

-  -  DF(Cov(T,T) ,rl,2,r2)  -  2*DF(Cov(T,T) ,rl,2)/r2 

-  -  DF(Cov(T,T) , rl , 2 , r2 , 2)  -  2*DF(Cov(T,T) ,rl,2,r2)/r2 

-  +  DF(Cov(T,T) , rl,2,r2,2) 

-  -  DF(Cov(T,T) , rl , 2 ,p2)/(G*R) 

-  -  DF(Cov(T,T) , rl,2,12)/(G*R*cos(p2)) 

-  -  DF(Cov(T,T) ,rl,2,r2)*DIRAC(r2  -  R')/(4*Pi*k) 


2 

Cov(6 , 1)  -  -  DF(Cov(T,T) ,pl)/(G  *R) 

Cov(6 , 2)  -  DF(Cov(T,T) ,pl,r2)/(G*R) 

Cov(6 , 3)  -  (DF(Cov(T,T) ,pl,r2)  +  2*DF(Cov(T,T) ,pl)/r2)/(G*R) 

Cov(6 ,4)  -  (DF(Cov(T,T) ,pl,r2,2)  +  2*DF(Cov(T,T) ,pl,r2)/r2)/(G*R) 
Cov(6 , 5)  -  -  DF(Cov(T,T)1pl,r2,2)/(G*R) 

2  2 

Cov(6.6)  -  DF(Cov(T, T) , pi , p2)/(G  *R  ) 


2  2 

Cov(6 , 7)  -  DF(Cov(T,T) ,pl , 12)/(G  *R  *cos(p2)) 

Cov(6 , 8)  -  DF(Cov(T,T) ,pl,r2)*DIRAC(r2  -  R' )/(4*Pi*k*G*R) 


2 

Cov(7 , 1)  -  -  DF(Cov(T,T),ll)/(G  *R*cos(pl)) 

Cov( 7,2)  -  DF(Cov(T,T) , ll,r2)/(G*R*cos(pl)) 

Cov( 7,3)  -  (DF(Cov(T,T) , 11 , r2)  +  2*DF(Cov(T,T) , ll)/r2)/(G*R*cos(pl) ) 
Cov(7,4)  -  (DF(Cov(T,T) , 11 ,r2 , 2)  +  2*DF(Cov(T,T) , 11 , r2)/r2)/(G*R*cos (pi) ) 
Cov( 7,5)  -  -  DF(Cov(T,T) , 11 , r2 , 2 )/(G*R*cos(pl) ) 

2  2 

Cov( 7,6)  -  DF(Cov(T,T),ll,p2)/(G  *R  *cos(pl)) 


2  2 

Cov( 7,7)  -  DF(Cov(T,T) ,11,12)/(G  *R  *cos(pl)*cos(p2) ) 

Cov( 7,8)  -  DF(Cov(T,T) , 11 ,r2)*DIRAC(r2  -  R ' ) / ( 4*Pi*k*G*R*cos (pi)) 
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CovC 8,1)  -  -  DF(Cov(T,T) ,rl)*DIRAC(rl  -  R' )/(4*Pi*k*G) 


Cov(8 , 2)  -  DF(Cov(T , T) , r 1 , r2)*DIRAC(rl  -  R')/(4*Pi*k) 

Cov( 8 , 3 )  -  DIRAC(rl  -  R' )*(DF(Cov(T,T) ,rl,r2)  +  2*DF(Cov(T,T) ,rl)/r2)/(4*Pi*k) 

Cov(8 ,4)  -  DIRAC(rl  -  R' )*(DF(Cov(T ,T) , rl , r2 , 2)  +  2*DF(Cov(T,T) , rl , r2)/r2)/ 
(4*Pi*k) 

Cov(8 , 5)  -  -  DF(Cov(T , T) , rl , r2 , 2)*DIRAC(rl  -  R’)/(4*Pi*k) 

CoW8,6)  -  DF(Cov(T , T) , rl , p2)*DIRAC(rl  -  R' )/(4*Pi*k*G*R) 

Cov(8,7)  -  DF(Cov(T , T) , rl , 12)*DIRAC(rl  -  R' )/(4*Pi*k*G*R*cos(p2)) 

2  2 

Cov(8 ,8)  -  DF(Cov(T,T) , rl , r2)*DIRAC(rl  -  R')*DIRAC(r2  -  R')/(16*Pi  *k  ) 


To  get  the  derivatives  r  '  r’:  e  autocovariance  function  Cov(T,T)  of  the 
disturbing  potential,  first  the  derivatives  of  the  expressions  B,  C 
and  Psi  have  to  be  com, ated.  This  was  also  done  by  means  of  REDUCE: 


2  2  4 

Derivatives  of  B  -  rl  *  r2  -  (R-D) 


2 

DF(B,r2)  -  2*rl  *r2 
2 

DF(B, r2 , 2)  -  2*rl 


2 

DF( B , rl)  -  2*rl*r2 
DF(B , rl , r2)  -  4*rl*r2 
DF(B , rl , r2 , 2)  -  4*rl 


2 

DF( B , r 1 , 2)  -  2*r2 
DF(B, rl , 2 , r2)  -  4*r2 
DF(B , r 1 , 2 , r2 , 2)  -  4 
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Derivatives  of  Psi  : 


cos(Psi)  -  sin(pl)  *  sin(p2)  +  cos(pl)  *  cos(p2)  *  cos(12-ll) 

DF(Psi ,p2)  -  (  -  sin(pl)*cos(p2)  +  sin(p2)*cos(ll  -  12)*cos(pl) )/sin(Psi) 
DF(Psi , 12)  -  -  sin(ll  -  12)*cos(pl)*cos(p2)/sin(Psi) 

DF(Psi,pl)  -  (sin(pl)*cos(ll  -  12)*cos(p2)  -  sin(p2)*cos(pl))/sin(Psi) 

2  2 

DF(Psi,pl,p2)  -  (cos(Psi)*(sin(pl)  *cos(ll  -  12)*cos(p2) 

2 

-  sin(pl)*sin(p2)*cos(ll  -  12)  *cos(pl)*cos(p2) 

-  sin(pl)*sin(p2)*cos(pl)*cos(p2) 

2  2 
+  sin(p2)  *cos(ll  -  12)*cos(pl)  ) 

2 

-  sin(Psi)  *(sin(pl)*sin(p2)*cos(ll  -  12) 

3 

+  cos(pl)*cos(p2)))/sin(Psi) 

DF(Psi,pl,12)  -  (cos(Psi)*sin(ll  -  12)*cos(pl)*cos(p2) 

*(sin(pl)*cos(ll  -  12)*cos(p2)  -  sin(p2)*cos(pl) ) 

2  3 

+  sin(Psi)  *sin(ll  -  12)*sin(pl)*cos(p2) )/sin(Psi) 

DF(Psi.ll)  -  sin(ll  -  12)*cos(pl)*cos(p2)/sin(Psi) 

DF(Psi , 11 ,p2)  -  (cos(Psi)*sin(ll  -  12)*cos (pl)*cos(p2) 

*(sin(pl)*cos(p2)  -  sin(p2)*cos(ll  -  12)*cos(pl)) 

2  3 

-  sin(Psi)  *sin(ll  -  12)*sin(p2)*cos(pl))/sin(Psi) 

2  2  2 

DF(Psi , 11 , 12)  -  (cos(Psi)*sin(ll  -  12)  *cos(pl)  *cos(p2) 

2  3 

-  sin(Psi)  *cos(ll  -  12)*cos(pl)*cos(p2))/sin(Psi) 
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2  2  4  2 

Derivatives  of  C  •  rl  *  r2  +  (R-D)  -  2*rl*r2*(R-D)  *cos(Psi)  : 


2 

DF(C,p2)  -  2*DF(Psi ,p2)*sin(Psi)*rl*r2*(R-D) 

2 

DF(C,12)  -  2*DF(Psi , 12)*sin(Psi)*rl*r2*(R-D) 

2  2 

DF(C,r2)  -  -  2*cos(Psi)*rl*(R-D)  +  2*rl  *r2 

2 

DF(C , r2 , 2)  -  2*rl 


2 

DF(C,pl)  -  2*DF(Psi ,pl)*sin(Psl)*rl*r2*(R-D) 

2 

DF(C,pl,p2)  -  2*cos(Psi)*DF(Psi , pl)*DF(Psl,p2)*rl*r2*(R-D) 

2 

+  2*DF(Psi , pi ,p2)*sin(Psi)*rl*r2*(R-D) 

2 

DF(C , pi , 12)  -  2*cos(Psi)*DF(Psi ,pl)*DF(Psi , 12)*rl*r2*(R-D) 

2 

+  2*DF(Psi ,pl , 12)*sin(Psi)*rl*r2*(R-D) 

2 

DF(C,pl,r2)  -  2*DF(Psi,pl)*sin(Psi)*rl*(R-D) 

DF(C,pl,r2,2)  -  0 


2 

DF(C, 11)  -  2*DF(Psi , ll)*sin(Psi)*rl*r2*(R-D) 

2 

DF(C, 11 ,p2)  -  2*cos(Psi)*DF(Psi , ll)*DF(Psi ,p2)*rl*r2*(R-D) 

2 

+  2*DF(Psl , 11 ,p2)*sln(Psl)*rl*r2*(R-D) 

2 

DF(C ,11,12)  -  2*cos(Psi)*DF(Psi , ll)*DF(Psl , 12)*rl*r2*(R-D) 

2 

+  2*DF(Psi , 11 , 12)*sln(Psi)*rl*r2*(R-D) 
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2 

DF(C,ll,r2)  -  2*DF(Psi , ll)*sin(Psi)*rl*(R-D) 
DF(C , 11 , r2 , 2)  -  0 


2  2 

DF(C , rl)  -  -  2*cos(Psi)*r2*(R-D)  +  2*rl*r2 


2 

DF(C,rl,p2)  -  2*DF(Psi ,p2)*sin(Psi)*r2*(R-D) 

2 

DF(C , rl , 12)  -  2*DF(Psi , 12)*sin(Psi)*r2*(R-D) 

2 

DF(C,rl,r2)  -  -  2*cos(Psi)*(R-D)  +  4*rl*r2 

DF(C , rl , r2 , 2)  -  4*rl 


DF(C,rl,2)  -  2*r2 
DF(C,rl,2,p2)  -  0 
DF(C , rl , 2 , 12)  -  0 
DF(C,rl,2,r2)  -  4*r2 
DF(C,rlI2,r2,2)  -  4 


Derivatives  of  Cov(T.T) 


5/2 

DF(Cov(T,T) ,p2)  -  -  3*DF(C,p2)*A*B*Co(T)/(2*C  ) 

5/2 

DF(Cov(T,T) ,12)  -  -  3*DF(C, 12)*A*B*Co(T)/(2*C  ) 

5/2 

DF(Cov(T,T) ,r2)  -  A*Co(T)*(2*DF(B,r2)*C  -  3*DF(C,r2)*B)/(2*C  ) 

2 

DF(Cov(T,T) , r2 , 2)  -  A*Co(T)*(4*DF(B,r2,2)*C  -  12*DF(B,r2)*DF(C,r2)*C 

2  7/2 

-  6*DF(C,r2,2)*B*C  +  15*DF(C,r2)  *B)/(4*C  ) 
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5/2 

DF(Cov(T,T) ,pl)  -  -  3*DF(C,pl)*A*B*Co(T)/(2*C  ) 


)p 

DF(Cov(T,T) ,pl,p2)  -  3*A*B*Co(T)*(  -  2*DF(C ,pl , p2 )*C  +  5*DF(C , pl)*DF(C , p2) )/(4*C 

1/2 

DF(Cov(T,T) , pi , 12)  -  3*A*B*Co(T)*(  -  2*DF(C , pi , 12)*C  +  5*DF(C, pl)*DF(C, 12) )/(4*C 
DF(Cov(T,T) ,pl,r2)  -  3*A*Co(T)*(  -  2*DF(B,r2)*DF(C,pl)*C  -  2*DF(C,pl , r2)*B*C 

7/2 

+  5*DF(C,pl)*DF(C,r2)*B)/(4*C  ) 


DF(Cov(T,T) ,pl,r2,2)  -  3*A*Co(T) 

2 

*(  -  4*DF(B,r2,2)*DF(C,pl)*C 


2 

8*DF(B,r2)*DF(C,pl,r2)*C 


2 

+  20*DF(B,r2)*DF(C,pl)*DF(C,r2)*C  -  4*DF(C ,pl , r2 , 2 )*B*C 


+  20*DF(C,pl,r2)*DF(C,r2)*B*C 

2 

+  10*DF(C,pl)*DF(C,r2,2)*B*C  -  35*DF(C,pl)*DF(C, r2)  *B) 
9/2 

/(8*C  ) 


5/2 

DF(Cov(T,T) ,11)  -  -  3*DF(C , ll)*A*B*Co(T)/(2*C  ) 

DF(Cov(T,T) , ll,p2)  -  3*A*B*Co(T) 

7/2 

*(  -  2*DF(C,ll,p2)*C  +  5*DF(C , 11)*DF(C ,p2) )/ (4*C  ) 

DF(Cov(T,T) ,11,12)  -  3*A*B*Co(T) 

7/2 

*(  -  2*DF(C,11,12)*C  +  5*DF(C,11)*DF(C,12))/(4*C  ) 

DF(Cov(T,T),ll,r2)  -  3*A*Co(T)*(  -  2*DF(B,r2)*DF(C, 11)*C  ^ 

-  2*DF(C,ll,r2)*B*C  +  5*DF(C,ll)*DF(C,r2)*B)/(4*C  ) 


DF(Cov(T , T) , 11 , r2 , 2)  -  3*A*Co(T) 

2  2 
*(  -  4*DF(B,r2,2)*DF(C,ll)*C  -  8*DF(B,r2)*DF(C,ll,r2)*C 

2 

-  4*DF(C, 11, r2, 2)*B*C  +  20*DF(C,ll,r2)*DF(C,r2)*B*C 

2 

+  10*DF(C, ll)*DF(C,r2,2)*B*C  -  35*DF(C, 11)*DF(C, r2)  *B 

9/2 

+  20*DF(B,r2)*DF(C,ll)*DF(C,r2)*C)/(8*C  ) 
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5/2 

DF(Cov(T,T) , rl)  -  A*Co(T)*(2*DF(B,rl)*C  -  3*DF(C,rl)*B)/(2*C  ) 


7/2 


DF(Cov(T,T) ,rl,p2)  -  3*A*Co(T)*<  -  2*DF(B,rl)*DF(C,p2)*C 

-  2*DF(C,rl,p2)*B*C  +  5*DF(C,rl)*DF(C,p2)*B)/(4*C  ) 

DF(Cov(T,T) ,rl,12)  -  3*A*Co(T)*(  -  2*DF(B,rl)*DF(C,12)*C 

7/2 

-  2*DF(C, rl , 12)*B*C  +  5*DF(C, rl)*DF(C, 12)*B)/(4*C  ) 

DF(Cov(T,T) ,rl,r2)  -  A*Co(T) 

2 

*(4*DF(B,rl,r2)*C  -  6*DF(B,rl)*DF(C,r2)*C 

-  6*DF(B,r2)*DF(C,rl)*C  -  6*DF(C,rl,r2)*B*C 

7/2 

+  15*DF(C, rl)*DF(C, r2)*B)/(4*C  ) 

DF(Cov(T,T) , rl , r2 , 2)  -  A*Co(T) 

3  2 

*(8*DF(B,rl,r2,2)*C  -  24*DF(B,rl,r2)*DF(C,r2)*C 

2  2 

-  12*DF(B,rl)*DF(C,r2 , 2)*C  +  30*DF(B,rl)*DF(C, r2)  *C 

2  2 

-  12*DF(B,r2,2)*DF(C,rl)*C  -  24*DF(B,r2)*DF(C,rl,r2)*C 


+  60*DF(B,r2)*DF(C,rl)*DF(C,r2)*C  -  12*DF(C,rl,r2,2)*B*C 
+  60*DF(C,rl,r2)*DF(C,r2)*B*C 

2 

+  30*DF(C,rl)*DF(C,r2 ,2)*B*C  -  105*DF(C,rl)*DF(C,r2)  *B) 
9/2 

/ (8*C  ) 


2 

DF(Cov(T,T) , rl , 2)  -  A*Co(T)*(4*DF(B, rl , 2)*C  -  12*DF(B,rl)*DF(C,rl)*C 

2  7/2 

-  6*DF( C , r 1 , 2 ) *B*C  +  15*DF(C,rl)  *B)/(4*C  ) 
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DF(Cov(T,T) ,rl,2,p2)  -  3*A*Co(T)  ^ 

*(  -  4*DF(B,rl,2)*DF(C,p2)*C  -  8*DF(B, rl)*DF(C ,rl ,p2)*C 

+  20*DF(B, rl)*DF(C, rl)*DF(C,p2)*C 

2 

+  20*DF(C, rl ,p2)*DF(C, rl)*B*C  -  4*DF(C , rl , 2 ,p2)*B*C 

2 

+  10*DF(C,rl,2)*DF(C,p2)*B*C  -  35*DF(C,rl)  *DF(C,p2)*B) 
9/2 

/(8* C  ) 

DF(Cov(T,T) ,rl,2,12)  -  3*A*Co(T) 

2  2 
*(  -  4*DF(B , rl , 2 )*DF(C , 12)*C  -  8*DF(B , rl)*DF(C , rl , 12)*C 

+  20*DF(B, rl)*DF(C, rl)*DF(C , 12)*C 

2 

+  20*DF(C,rl,12)*DF(C,rl)*B*C  -  4*DF(C,rl,2,12)*B*C 

2 

+  10*DF(C , rl , 2)*DF(C , 12)*B*C  -  35*DF(C,rl)  *DF(C,12)*B) 
9/2 

/(8 *C  ) 

DF(Cov(T,T) .rl,2,r2)  -  A*Co(T) 

2  3 

*(  -  24*DF(B,rl,r2)*DF(C,rl)*C  +  8*DF(B,rl,2,r2)*C 

2  2 

-  12*DF(B, rl, 2)*DF(C, r2)*C  -  24*DF(B,rl)*DF(C,rl,r2)*C 

+  60*DF(B , rl)*DF(C , rl)*DF(C , r2 )*C 

2  2 

-  12*DF(B,r2)*DF(C,  rl,2)*C  +  30*DF(B,r2)*DF(C, rl)  *C 

2 

+  60*DF(C, rl , r2)*DF(C, rl)*B*C  -  12*DF(C,rl,2,r2)*B*C 

2 

+  30*DF(C,rl,2)*DF(C,r2)*B*C  -  105*DF(C,rl)  *DF(C,r2)*B) 
9/2 

/(8*C  ) 
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DF(Cov(T,T),rl.2,r2.2)  -  A*Co<T) 


*(  -  48*DF(B,rl,r2,2)*DF<C,rl)*C 

3 

-  96*DF(B,rl,r2)*DF(C,rl,r2)*C 

2 

+  240*DF(B,rl,r2)*DF(C,rl)*DF(C,r2)*C 

4 

+  16*DF(B,rl, 2 ,r2 ,2)*C 

3 

-  48*DF(B,rl,2,r2)*DF(C,r2)*C 

3 

-  24*DF(B,rl,2)*DF(C,r2,2)*C 

2  2 

+  60*DF(B,rl,2)*DF(C,r2)  *C 

3 

-  48*DF(B,rl)*DF(C,rl,r2,2)*C 

2 

+  240*DF(B,rl)*DF(C,rl,r2)*DF(C,r2)*C 

2 

+  120*DF(B,rl)*DF(C,rl)*DF(C,r2 ,2)*C 

2 

-  420*DF(B,rl)*DF(C,rl)*DF(C,r2)  *C 

3 

-  24*DF(B,r2 , 2)*DF(C, rl , 2)*C 

2  2 

+  60*DF(B,r2,2)*DF(C,rl)  *C 

2 

+  240*DF ( B , r 2 ) *DF ( C , rl , r2)*DF(C , rl)*C 

3 

-  48*DF(B,r2)*DF(C,rl,2,r2)*C 

2 

+  120*DF(B,r2)*DF(C,rl,2)*DF(C,r2)*C 

2 

-  420*DF(B,r2)*DF(C,rl)  *DF(C,r2)*C 

2 

+  120*DF(C, rl , r2 , 2)*DF(C , rl)*B*C 
2  2 

+  120*DF(C,rl,r2)  *B*C 

-  840*DF(C,rl,r2)*DF(C,rl)*DF(C,r2)*B*C 

3 

-  24*DF(C, rl , 2 ,r2 , 2)*B*C 

2 

+  120*DF(C , rl , 2 , r2)*DF(C , r2)*B*C 

2 

+  60*DF(C , rl , 2)*DF(C , r2 , 2)*B*C 

2 

-  210*DF(C, rl , 2)*DF(C , r2)  *B*C 

2 

-  210*DF(C,rl)  *DF(C,r2,2)*B*C 

2  2 

+  945*DF(C,rl)  *DF(C,r2)  *B) 


11/2 
/(16*C  ) 
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Asymptotic  relations  (flat  earth  model) 

For  local  applications  the  asymptotic  form  (6-21)  is  used  for  the  autocovariance 
function  of  the  disturbing  potential 

2  A 

Cov(T,T)(rl,r2,Psi)  =  4  D  .  Co(T)  , 

3/2 

B 

with 

A  -  2D  +  zl  +  z2 

2  2 
B  -  S  +  A 

2  2  2 
S  -  U  +  V 

U  -  x2  -  xl 

V  -  y2  -  yl  . 

For  the  functionals  of  the  gravity  field  in  P  we  get  the  following  asymptotic  forms 

1)  geoidal  height  N1  -  Tl/G 

2)  gravity  disturbance  d  gl  -  -  DF(Tl.zl) 

3)  gravity  anomaly  D  gl  -  -  DF(Tl.zl) 

4)  radial  gravity  gradient  -  -  DF(Tl,zl,2) 

5)  second  radial  derivative  of  T  -  +  DF(Tl,zl,2) 

6)  N-S -component  of  the 

deflection  of  the  vertical  Xil  -  -  DF(Tl,xl)/G 

/)  K-W- component  of  the 

deflection  of  the  vertical  Etal  -  -  DF(Tl,yl)/G 

8)  density  contrast  d  rho  -  -  DF(T1 , zl)*DIRAC(zl)/(4*Pi*k) 

To  get  the  functionals  in  Q  substitute  Tl,xl,yl,zl  by  T2,x2,y2,z2. 

Noteworthy  is  that  the  covariances  of  the  gravity  disturbance  and  the  gravity 
anomaly  become  identical. 
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The  covariances  between  the  functionals  are 

2 

Cov(l , 1)  -  +  Cov(T.T)  /  G 
Cov(l,2)  -  -  DF(Cov(T,T),z2)  /  G 
Cov(l,3)  -  -  DF(Cov(T,T),z2)  /  G 
Cov(l ,4)  -  -  DF(Cov(T,T) ,z2,2)  /  G 
Cov(l , 5)  -  +  DF(Cov(T,T) ,z2,2)  /  G 

2 

Cov(l , 6)  -  -  DF(Cov(T,T) ,x2)  /  G 

2 

Cov(l,7)  -  -  DF(Cov(T,T),y2)  /  G 

Cov(l , 8)  -  -  DF(Cov(T,T) , z2)  *  DIRAC(z2)  /  (4*Pi*k*G) 

Cov(2 , 1)  -  -  DF(Cov(T,T),zl)  /  G 
Cov(2 , 2)  -  +  DF(Cov(T,T) ,zl,z2) 

Cov(2,3)  -  +  DF(Cov(T,T) ,zl,z2) 

Cov(2 , 4)  -  +  DF(Cov(T,T),zl,z2,2) 

Cov(2 , 5)  -  -  DF(Cov(T,T) ,zl,z2,2) 

Cov(2 ,6)  -  +  DF(Cov(T,T) ,x2 ,zl)  /  G 
Cov(2 , 7)  -  +  DF(Cov(T,T) ,y2 , zl)  /  G 

Cov(2 , 8)  -  +  DF(Cov(T,T) ,zl,z2)  *  DIRAC(z2)  /  (4*Pi*k) 

Cov(3,l)  -  -  DF(Cov(T.T) ,zl)  /  G 
Cov(3 , 2)  -  +  DF(Cov(T,T),zl,z2) 

Cov(3 , 3)  -  +  DF(Cov(T,T) ,zl,z2) 

Cov(3 ,4)  -  +  DF(Cov(T,T) , zl , z2 , 2) 

Cov(3 , 5)  -  -  DF(Cov(T,T) ,zl,z2,2) 

Cov(3,6)  -  +  DF(Cov(T,T) ,x2 , zl)  /  G 
Cov(3 , 7)  -  +  DF(Cov(T,T) ,y2 ,zl)  /  G 

Cov(3,8)  -  +  DF(Cov(T,T),2l,z2)  *  DIRAC(z2)  /  (4*Pi*k) 


134 


Cov(4,l)  - 
Cov(4,2)  - 
Cov(4,3)  - 
Cov(4,4)  — 
Cov(4,5)  - 
Cov(4,6)  - 
Cov(4,7)  - 
Cov(4,8)  - 

Cov( 5,1)  - 
Cov( 5 , 2 )  - 
Cov(5 , 3)  - 
Cov(5,4)  - 
Cov(5,5)  - 
Cov(5,6)  - 
Cov( 5 , 7 )  - 
Cov( 5,8)  - 

Cov(6,l)  - 
Cov(6,2)  - 
Cov(6,3)  - 
Cov(6,4)  - 
Cov(6,5)  - 

Cov(6,6)  ~ 

Cov(6 , 7 )  - 
Cov(6,8)  - 


-  DF(Cov(T , T) , zl , 2)  /  G 
+  DF ( Cov(T ,T) , zl , 2 , z2) 

+  DF(Cov(T ,T) ,zl,2,z2) 

+  DF(Cov(T, T) , zl , 2 , z2 , 2) 

-  DF(Cov(T.T) ,zl,2,z2,2) 

+  DF(Cov(T , T) , x2 , zl , 2)  /  G 
+  DF(Cov(T,T) , y2 , zl , 2)  /  G 

+  DF(Cov(T ,T) , zl , 2 , z2)  *  DIRAC(z2)  /  (4*Pi*k) 

+  DF(Cov(T , T) , zl , 2)  /  G 

-  DF(Cov(T ,T) , zl , 2 , z2) 

-  DF(Cov(T , T) , zl , 2 , z2) 

-  DF(Cov(T,T) ,zl,2,z2,2) 

+  DF(Cov(T ,T) , zl , 2 ,z2 , 2) 

-  DF(Cov(T,T) ,x2 ,zl , 2)  /  G 

-  DF(Cov(T,T) ,y2,zl,2)  /  G 

-  DF(Cov(T,T) ,zl,2,z2)  *  DIRAC(z2)  /  (4*Pi*k) 

2 

-  DF(Cov(T,T) ,xl)  /  G 

+  DF(Cov(T,T) ,xl,z2)  /  G 
+  DF(Cov(T,T) , xl,z2)  /  G 
+  DF(Cov(T,T),xl,z2,2)  /  G 

-  DF(Cov(T,T) ,xl,z2 ,2)  /  G 

2 

+  DF(Cov(T,T) ,xl,x2)  /  G 

2 

+  DF(Cov(T,T) ,xl,y2)  /  G 

+  DF(Cov(T,T),xl.z2)  *  DIRAC(z2)  /  (4*Pi*k*G) 
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Cov(7 , 1) 
Cov(7 ,2) 
Cov(7 , 3) 
Cov(7 ,4) 
Cov(7 , 5) 

Cov(7 ,6) 

Cov(7 ,7) 
Cov(7 , 8) 


2 

-  DF(Cov(T ,T) ,yl)  /  G 

+  DF(Cov(T,T) ,yl,z2)  /  G 
+  DF(Cov(T,T) ,yl,z2)  /  G 
+  DF(Cov(T,T) ,yl , z2 , 2)  /  G 

-  DF(Cov(T,T) ,yl , z2 , 2)  /  G 

2 

+  DF(Cov(T,T) ,x2 ,yl)  /  G 

2 


+  DF(Cov(T,T),yl,y2)  /  G 

+  DF(Cov(T,T),yl,z2)  *  DIRAC(z2)  /  (4*Pi*k*G) 


Cov(8 , 1) 
Cov (8,2) 
Cov(8 ,3) 
Cov(8 ,4) 
Cov(8 , 5) 
Cov(8,6) 
Cov(8 , 7) 
Cov(8 , 8) 


-  DF(Cov(T,T) ,zl)  *  DIRAC(zl)  /  (4*Pi*k*G) 

+  DF(Cov(T,T) ,zl,z2)  *  DIRAC(zl)  /  (4*Pl*k) 

+  DF(Cov(T,T) ,zl,z2)  *  DIRAC(zl)  /  (4*Pi*k) 

+  DF(Cov(T,T) ,zl,z2,2)  *  DIRAC(zl)  /  (4*Pi*k) 

-  DF(Cov(T,T) ,zl,z2,2)  *  DIRAC(zl)  /  (4*Pl*k) 

+  DF(Cov(T,T) ,x2 ,zl)  *  DIRAC(zl)  /  (4*Pi*k*G) 

+  DF(Cov(T,T) ,y2,zl)  *  DIRAC(zl)  /  (4*Pl*k*G) 

2  2 

+  DF(Cov(T,T) , zl , z2)  *  DIRAC(zl)  *  DIRAC(z2)  /  (16*P1  *k  ) 


In  this  case  the  derivatives  of  the  autocovariance  function  Cov(T,T)  of  the 
disturbing  potential  can  be  given  directly  as 


2  5/2 

DF(Cov(T,T) ,x2)  -  -  12*D  *  Co(T)  *  A  *  U  /  B 

2  5/2 

DF(Cov(T,T) , y 2 )  -  -  12*D  *  Co(T)  *  A  *  V  /  B 

2  2  5/2 

DF(Cov(T,T) ,z2)  -  4*D  *  Co(T)  *  (  -  3*A  +  B)  /  B 

2  2  7/2 

DF(Cov(T,T),z2,2)  -  12*D  *  Co(t)  *  A  *  (5*A  -  3*B)  /  B 
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2  5/2 

DF(Cov(T,T) ,xl)  -  12*D  *  Co(T)  *  A  *  U  /  B 


2  2  7/2 

DF(Cov(T,T),xl,x2)  -  12*D  *  Co(T)  *  A  *  (B  -  5*U  )  /  B 

2  7/2 

DF(Cov(T,T) ,xl,y2)  -  -  60*D  *  Co(T)  *A*U*V/B 

2  2  7/2 

DF(Cov(T,T) ,xl,z2)  -  12*D  *  Co(T)  *  U  *  (  -  5*A  +  B)  /  B 

2  2 

DF(Cov(T,T) ,xl,z2,2)  -  60*D  *  Co(T)  *  A  *  U  *  (7*A  -  3*B)  / 


2  5/2 

DF(Cov(T,T) ,yl)  -  12*D  *  Co(T)  *  A  *  V  /  B 

2  7/2 

DF(Cov(T,T) ,yl,x2)  -  -  60*D  *  Co(T)  *A*U*V/B 

2  2  7/2 

DF(Cov(T,T) ,yl,y2)  -  12*D  *  Co(T)  *  A  *  (B  -  5*V  )  /  B 

2  2  7/2 

DF(Cov(T,T) ,yl ,z2)  -  12*D  *  Co(T)  *  V  *  (  -  5*A  +  B)  /  B 

2  2 

DF(Cov(T,T) ,yl,z2,2)  -  60*D  *  Co(T)  *  A  *  V  *  (7*A  -  3*B)  / 


2  2  5/2 

DF(Cov(T,T) ,zl)  -  4*D  *  Co(T)  *  (  -  3*A  +  B)  /  B 

2  2  7/2 

DF(Cov(T,T) ,zl,x2)  -  12*D  *  Co(T)  *  U  *  (5*A  -  B)  /  B 

2  2  7/2 

DF(Cov(T,T) ,zl,y2)  -  12*D  *  Co(T)  *  V  *  (5*A  -  B)  /  B 

2  2  7/2 

DF(Cov(T,T) ,zl,z2)  -  12*D  *  Co(T)  *  A  *  (5*A  -  3*B)  /  B 

2  4  2 

DF(Cov(T,T) ,zl,z2,2)  -  12*D  *  Co(T)  *  (  -  35*A  +  30*A  *B  - 


9/2 

B 


9/2 

B 


2 

3*B  )  /  B 


2  2  7/2 

DF(Cov(T,T) , zl , 2)  -  12*D  *  Co(T)  *  A  *  (5*A  -  3*B)  /  B 


2 

2 

9/2 

DF(Cov(T,T) , zl,2,x2)  -  60*D 

*  Co(T)  *  A  *  U  *  ( 

-  7*A  +  3*B) 

/  B 

2 

2 

9/2 

DF(Cov(T,T) ,zl,2,y2)  -  60*D 

*  Co(T)  *  A  *  V  *  ( 

-  7*A  +  3*B) 

/  B 

2  4  2  2  9/2 

DF(Cov(T,T) , zl , 2 , z2)  -  12*D  *  Co(T)  *  (  -  35*A  +  30*A  *B  -  3*B  )  /  B 

2  4  2  2  11/2 

DF(Cov(T,T) , zl , 2 , z2 , 2)  -  60*D  *  Co(T)  *  A  *  (63*A  -  70*A  *B  +  15*B  )  /  B 
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